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Introduction  and  Executive  Summary 


The  trend  towards  very  large  area  ratio  nozzles,  which  result  in  performance  gains  for 
space  propulsion  applications,  has  increased  the  need  for  detailed  knowledge  of  the  momentum 
losses  due  to  nozzle  viscous  effects  (i.e.,  boundary  layer).  These  losses  degrade  overall  system 
performance,  such  as  increasing  system  weight,  decreasing  useful  payload  weight,  and/or 
decreasing  effective  system  range.  Another  important  factor  in  the  designing  of  propulsive 
nozzles  is  the  detailed  knowledge  of  heat  transfer  at  the  wall  for  regeneratively  cooled  walls 
and/or  material  performance. 

Because  of  the  importance  to  rocket  propulsion,  the  Astronautics  Laboratory 
(AL/AFSC)  has  sponsored  the  Boundary  Layer  Study  Contract  to  improve  the  understanding 
and  computational  predictive  capabilities  for  boundary  layers  in  rocket  nozzles  with  very  high 
area  ratios. 

The  contractual  effort  was  broken  down  into  4  work  tasks  with  2  deliverable  computer 
programs  plus  documentation  and  two  technical  reports.  The  basic  thrust  of  the  effort  was  to 
develop  an  increased  analytical  capability  to  predict  the  performance  loss  of  thick  boundary 
layers  in  high  area  ratio  propulsive  nozzles.  This  objective  has  been  successfully  met. 

A  description  of  the  four  work  tasks  under  this  contract,  along  with  the 
accomplishments  made  during  each  task,  are: 

1)  Extension  of  Current  Boundary  Layer  Technology: 

o  Extended  the  Boundary  Layer  capability  of  the  JANNAF  TDK/BLM 
code  to  predict  thick  boundary  layer  losses  including  the  effects  of 
longitudinal  curvature 

o  Established  the  deficiency  of  the  current  JANNAF  method  of  computing 
the  boundary  layer  thrust  loss  for  thick  boundary  layers 
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2)  Thick  Boundary  Layer  Assessment 


o  Surveyed  existing  codes  (both  full  and  parabolized  Navier-Stokes 
solvers)  for  solutions  to  the  nozzle  flow  problem 

o  Assessed  the  effects  of  various  assumptions  on  computed  boundary  layer 
parameters 

o  Established  the  course  of  action  to  be  taken  for  the  development  of  the 
new  cede  developed  in  the  Phase  3  effort 

3)  New  Analytical  Methods 

o  Produced  a  new  code  for  predicting  the  performance  (including  very 
thick  boundary  layers)  of  propulsive  nozzles.  The  new  code  is  called 
VIPER  (Viscous  Performance  Evaluation  Routine) 

o  Validated  the  VIPER  against  data  and  the  existing  JANNAF  Standard 
Code,  TDK/BLM. 

4)  Experimental  Validation  Test  Plan 

o  Assessed  the  state  of  the  art  of  diagnostics  techniques  for  boundary  layer 
measurements  in  nozzles  flows 

o  Visited  and  discussed  methods  and  techniques  with  sources  in  industry 
and  government 

o  Made  recommendations  for  a.  test  plan  to  validate  the  turbulence  models 
used  in  the  codes  developed  under  this  contract 

A  series  of  4  technical  reports  were  written  which  describe  in  detail  all  of  the  technical 
effort  under  this  contract.  These  reports  are: 
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Phase  1  Effort 


AFAL-TR-87-031,  "Two-Dimensional  Kinetics  (TDK)  Nozzle  Performance 
Program  -  Thick  Boundary  Layer  Version",  February  1988 

Phase  2  Effort 

AL-TR-90-041,  "Boundary  Layer  Study  Thick  Boundary  Layer  Assessment" 
Phase  3  Effort 

AL-TR-90-042,  "Viscous  Interaction  Performance  Evaluation  Routine  for 
Nozzle  Flows  with  Finite  Rate  Chemistry" 

Phase  4  Effort 

AL-TR-90-043,  "Boundary  Layer  Study  Experimental  Validation  Test  Plan" 


Conclusions 

o  The  more  general  form  of  the  JANNAF  Boundary  Layer  Loss  Equation 
should  be  adopted  as  the  JANNAF  standard 

o  That  longitudinal  curvature  is  not  an  important  effect  for  conventional 
propulsive  nozzles 

o  That  PNS  codes,  such  as  the  VIPER  code,  are  required  to  supply  the 
necessary  accuracy  to  resolve  the  boundary  layer  and  core  flow  for 
propulsive  nozzles  with  thick  boundary  layers 

o  The  VIPER  code  represents  a  significant  advancement  in  nozzle 

performance  prediction  capability 

o  That  direct  experimental  measurement  of  the  boundary  layer  loss 

parameters  are  required  to  validate  the  analytical  loss  models 
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Recommendations 


o  That  the  VIPER  code  be  considered  as  a  JANNAF  standard  for  high  area 
ratio  propulsive  nozzles 

o  That  the  VIPER  code  be  extended  to  add: 
particulate  flows, 

real  nozzle  effects,  such  as  tangential  slot  injection, 
chamber/injector  models, 
better  shock  capturing  capability,  and 
a  more  robust  k-e  turbulence  model 

o  That  the  experimental  test  plan  recommended  in  the  Pha.se  4  work  be 
executed 
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2.0  Phase  I  Work  Effort  -  Improvements  to  Current  Technology 


The  purpose  of  the  Phase  I  work  effort  was  to  investigate  ways  in  which  the  current 
performance  prediction  methodologies  could  be  extended,  and  to  select  the  most  promising 
approaches  and  implement  them  in  the  current  performance  code(s).  The  objectives  of  the 

*  Phase  I  work  were  met.  Section  2.1  describes  the  investigation  of  methods  to  extend  current 
technology,  section  2.2  describes  the  actual  modifications  and  results,  and  section  2.3  describes 

x  the  results  and  conclusions. 

2.1  Investigation  Into  Extensions  of  Current  Technology 

At  the  start  of  the  Boundary  Layer  Study  Contract,  the  April  1985  Version  of  the  TDK1 

program  was  the  JANNAF  standard  code  for  computing  the  performance  of  liquid  propellant 

rocket  engine  nozzles.  This  code  represented  the  state  of  the  art  in  nozzle  performance 

prediction.  The  boundary  layer  module  incorporated  into  TDK  at  that  time  was  the  BLM 
12  3 

model  ’  ’ .  This  module  would  compute  the  boundary  layer  in  a  nozzle  assuming  a 
non-calorically  perfect  gas  for  either  laminar  or  turbulent  flow.  The  turbulence  model  used 
was  the  Cebeci-Smith  eddy  viscosity  model ,  which  also  included  the  effect  of  transverse 
curvature. 

Since  the  purpose  of  the  Phase  I  work  was  to  extend  current  technology  in  computing 
boundary  layer  losses,  it  was  important  to  understand  the  underlying  assumptions  which  are 
present  in  the  boundary  layer  equations  as  used  for  performance  calculations. 

* 

The  boundary  layer  equations  are  "deduced"  using  an  order  of  magnitude  argument  . 
Reference  4,  pp  39-44,  contains  an  excellent  presentation  of  the  TSL  (thin  shear  layer) 
equation  derivation.  Classical  derivations  of  the  TSL  equations,  e.g.,  Schlichting  ,  pp  128-131, 
Laudau  and  Lifshitz  ,  pp  145-147,  tend  to  ignore  the  fact  that  the  normal  pressure  gradient, 
d?  /  dy,  vanishes  not  because  the  shear  layer  is  thin  (8  /  i  «1),  but  because  they  have 

assumed  that  the  surface  is  relatively  flat,  i.e.,  that  the  radius  of  curvature,  R  ,  is  large,  or 

*  c 


There  are  more  mathematic  llv  rieorous  derivations  of  the  boundary  laver  eauations 

7  ... 

than  these,  for  example,  Van  Dyke  .  However,  the  choice  of  the  perturbation 
parameter  and  form  of  the  series  expansion  are  based  on  order  of  magnitude  analysis. 
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more  precisely  Rc  /  l  »1.  The  term  that  is  normally  neglected  is  the  "centrifugal  Torce" 
equation,  i.e., 

9P  _  pu2 

5y 

The  neglect  of  this  term  is  justified  when  R  It »  8  /  l  which  is  usually  the  case  in  propulsive 

c 

nozzles  (5  is  usually  very  small  when  i  is  small).  However,  when  the  shear  layer  thickens  to 

t 

the  point  that  8  /  i  =  0(1),  then  we  can  no  longer  assume  a  priori  that  the  normal  pressure 
gradient  has  a  negligible  effect  on  the  flow  in  thv,  shear  layer. 

There  are  two  effects  which  can  have  a  potentially  significant  impact  on  the  viscous 
thrust  loss  which  are  not  normally  accounted  for  in  classical  boundary  layer  theory: 

o  transverse  curvature 

o  longitudinal  curvature 

The  -effect  of  transverse  curvature  was  already  treated  in  the  BLM  module  and  is  not 
discussed  here.  In  the  following  sections,  the  effect  of  longitudinal  curvatui ,  will  be  examined 
in  light  of  its  impact  on  the  mean  flow  and  the  eddy-viscosity  models.  Also  not  included  in 
the  JANNAF  boundary  layer  codes  is  the  effects  of  wall  roughness.  These  effects  are  also 
discussed. 
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Longitudinal  Curvature  Effects 


The  effect  of  longitudinal  curvature  on  the  flow  results  (i.e.,  nozzle  performance) 
manifests  itself  in  several  ways,  the  most  important  are 

0  on  the  mean  flow  in  the  hour. '’ary  layer  (centrifugal  force  and  normal  stresses) 

0  in  the  turbulent  shear  stress  model  (eddy  viscosity  models) 

0  in  the  interaction  with  the  core  flow 

The  variation  in  mean  flow  due  to  the  centrifugal  force  balance  nas  been  estimated  by 

8 

Beddini  to  be  about  two  percent  in  momentum  thickness  on  the  J2  nozzle.  This  change  in 
thickness  translates  into  a  like  amount  in  thrust  loss.  Note  that  while  none  of  the  commonly 
used  boundary  layer  codes  are  accurate  to  within  5%,  (probably  not  even  10%),  it  is  important 
to  be  able  to  predict  the  trends  accurately  in  order  to  know  the  relative  amount  of  performance 
increase  (or  decrease)  if  the  nozzle  area  ratio  is  enlarged. 

4 

Cebeci  &  Bradshaw  discuss  the  mean  flow  terms  and  their  order  of  magnitude  relative 
to  the  other  boundary-layer  terms.  They  point  out  that  centrifugal  force  terms  are  generally 
intractable  and  normal  stresses  are  usually  negligible  in  weakly  turbulent  flows.  For  this 
reason  the  thin  boundary-layer  approach  is  usually  adopted  even  in  cases  where  it  isn't  strictly 
applicable. 

The  contribution  of  turbulence  to  the  normal  stresses  have  received  relatively  little 

9 

attention  in  the  literature.  Finley  analyzed  experimental  data  from  cooled-wall  and  adiabatic 
hypersonic  flow  nozzles.  He  attributed  the  large  variation  of  static  pressure  in  hypersonic 
nozzle  boundary-layer  flows  to  be  the  combined  effects  of  longitudinal  curvature  of  the  mean 
streamlines  and  the  increasing  importance,  as  the  Mach  number  rises,  of  the  Reynolds  stress 
contribution  to  the  total  normal  stress  perpendicular  to  the  wall.  He  concluded  that  at  its  peak 
value  the  Reynolds  stress  may  provide  a  normal  stress  contribution  equal  to  that  of  the  mean 
static  pressure. 

The  influence  of  longitudinal  curvature  on  the  mean  boundary-layer  flow  can  be 
modeled  in  several  ways.  The  first  way,  pointed  out  by  many  authors,  is  to  simply  include  the 
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centrifugal  force  term  in  the  y-momentum  equation,  i.e. 


l-P-s-lr 


where  R  is  the  longitudinal  radius  of  curvature  of  the  body. 


Two  rigorous  approaches  to  modeling  the  effects  of  longitudinal  curvature  (on  the 
mean  flow)  have  been  taken  in  the  literature.  These  are  the  singular  perturbation  analysis  of 
Van  Dyke  ’  and  the  metric  influence  of  curvature  method  of  Schultz-Grunow  and  Brener  . 

The  singular  perturbation  theory  is  too  complex  to  go  into  here.  Suffice  it  to  say  that 
starting  with  the  Navier— Stokes  equations  the  problem  is  solved  by  a  scheme  of  successive 
approximations  by  the  method  of  inner  and  outer  expansions.  Two  complementary  expansions 
are  constructed  simultaneously,  and  matched  in  their  overlap  region  of  common  validity.  To 
lowest  order  in  the  inner  region  the  Prandtl  boundary-layer  theory  is  recovered  while  the  ouier 
region  is  an  inviscid  flow.  Higher-order  solutions  in  each  region  give  the  perturbations  to  the 
boundary-layer  and  external  flows.  For  example,  the  boundary-layer  displacement  effect 
shows  up  in  both  the  inner  and  outer  regions.  The  perturbation  method  simplifies  the  problem 
because  viscous  effects  are  confined  to  the  inner  region  which  is  parabolic  and  hence  easier  to 
solve  than  the  full  (elliptic)  Navier— Stokes  equations. 

The  approach  of  Schultz-Grunow  and  Breuer11  is  considerably  different  from  the 
perturbation  approach.  They  adopt  a  curvilinear  system  of  coordinates  and  take  the  arc  length 
along  the  surface  as  the  x-coordinate.  This  brings  the  longitudinal  curvature  directly  into  the 
problem  since  any  differentiation  with  respect  to  x  carries  a  factor  of  the  ratio  of  curvature 
radii  (wall  position  to  actual  position).  Schultz-Grunow  and  Breuer  derive  the  boundary  layer 
equations  from  an  order  of  magnitude  analysis  performed  on  the  complete  Navier-Stokes 
equations.  They  obtain  the  following  set  of  equations: 


Continuity: 


l  +  jjjrO+Wv.O 


8 


x-momentum 


y-momentum 

where  k  =  1/R  .  These  equations  are  valid  for  concave  as  well  as  convex  surfaces  (k>0  for 

V 

convex  walls  and  k<0  for  concave  walls). 

The  authors  point  out  that  these  equations  contain  all  the  second-order  terms  of  the 

second-order  perturbation  theory  (Van  Dyke10),  namely  those  of  order  1  /{Re,  and  the  exact 
metric  influence  of  curvature.  An  additional  term  was  included  in  the  equations  as  well  (the 
only  justification  given  is  that  the  main  flow  is  then  an  analytic  solution  of  the  boundary-layer 
equations  [sic]).  Also  note  that  these  equations  reduce  to  the  classical  boundary-layer 
equations  when  k  =  0. 

12 

Cebeci,  Hirsch,  &  Whitelaw  analyzed  the  turbulent  boundary  layer  on  a  (convex) 

li 

longitudinally  curved  surface  using  the  mean  flow  equations  of  Schultz-Grunow  and  Breuer 

discussed  above.  In  their  treatment  they  replaced  the  laminar  viscosity  with  an  effective 

viscosity,  v  +  e  .  The  eddy  viscosity  (em)  was  specified  to  be  the  standard  Cebeci  and 

13  14 

Smith  eddy-viscosity  formulation  modified  by  Bradshaw's  correction  for  longitudinal 

curvature. 

13 

Cebeci  and  Smith  point  out  that  streamline  curvature  may  increase  or  decrease  the 

turbulent  mixing,  depending  on  the  degree  of  wall  curvature,  and  it  can  strongly  affect  the  skin 

friction  and  heat  transfer  coefficients.  They  also  included  the  streamwise  curvature  effect  into 

15 

the  eddy-viscosity  expression.  (This  will  be  discussed  later).  Eghlima  et.al.  used  a  similar 
approach. 

Longitudinal  curvature  has  an  effect  on  the  structure  of  turbulence  as  well.  In  a 
boundary  layer  on  a  convex  wall,  see  Figure  1,  the  centrifugal  forces  exert  a  small  stabilizing 
effect.  In  ontrast  with  that,  concave  walls  have  a  de-stabilizing  effect  due  to  instabilities 
known  as  Taylor-Gortler  vortices,  (see  Figure  2).  Physically,  the  fluid  element  in  a  curved 
stream  whose  angular  momentum  decreases  with  increasing  distance  from  the  center  of 


u  I + <I+ky>v  |  • * i + <1+ky>v  0 


+  vk[ 


du  ku 


dy '  1+kv 


] 


ltUJ  =  iiisx 
P  ay 
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Figure  1:  #low  by  Longitudinally  Curved  Walls 
(From  Reference  16) 


Figure 


Streamline  pattern  of  Taylor-Oortler  Vortices 
(From  Reference  16) 


curvature  is  unstable.  If  the  element  is  slightly  displaced  outward  from  the  center,  conserving 
its  ang.  lar  momentum  about  the  center,  it  will  be  moving  faster  than  its  surroundings  (which 
have  a  smaller  angular  momentum).  Therefore  the  radial  pressure  gradient  that  controls 
streamline  curvature  will  be  too  small  to  direct  the  displaced  element  along  a  streamline  of  the 
main  flow,  and  the  element  will  move  even  further  outward.  The  converse  argument  holds  for 
inward  displacement.  These  instabilities  affect  both  the  transition  to  turbulence  and  the  eddy 
viscosity. 

14 

Bradshaw  has  developed  an  expression  for  the  effect  of  longitudinal  curvature  which 

13 

can  be  incorporated  into  the  Cebeci-Smith  eddy  viscosity  expression  .  The  inner 
eddy-viscosity  expression  is  multiplied  by  S2  which  was  obtained  by  Bradshaw  using  an 
analogy  between  streamline  curvature  and  buoyancy  in  turbulent  shear  flow.  The  expression 
for  S  is  given  by 

S  =  1/(1  +  /3R.) 

Rj  =  2u/(Rc  du/<5y) 


where  Ri  is  analogous  to  the  Richardson  number.  The  parameter  /3  is  equal  to  7  for  convex 

7 

surfaces  and  4  for  concave  surface  according  to  meteorological  data  and  by  Bradshaw’s 
analogy  between  streamline  curvature  and  buoyancy.  According  to  Bradshaw,  the  effects  of 
curvature  on  the  mixing  length  or  eddy  viscosity  are  appreciable  if  the  ratio  of  boundary-layer 
thickness  to  radius  of  curvature,  8/R  ,  exceeds  roughly  1/300. 

V 


The  interaction  of  the  core  flow  with  the  shear  layer  manifests  itself  in  the  application 
of  the  boundary  conditions  for  both  the  shear  layer  flow  and  the  core  flow.  For  the  TSL 
approximation  it  is  proper  to  view  the  system  as  a  singular  perturbation  problem  with  the  inner 
flow  (i.e.,  the  shear  layer)  boundary  condition  being  applied  at  infinity  (not  8),  and  the  outer 
flow  (the  core  flow)  boundary  condition  applied  at  the  wall.  The  matching  condition  is: 

Inner  Limit  of  the  Outer  Expansion  (ILOE)  = 

Outer  Limit  of  the  Inner  Expansion  (OLIE) 


11 


The  two  flows  are  then  connected  iteratively  by  displacing  the  wall  by  8*  and  recomputing  the 
outer  flow. 


The  same  procedure  of  matching  the  two  flows  (i.e.,  ILOE=OLIE)  is  used  for  higher 

7 

order  boundary  layer  theory  as  discussed  earlier  (see  Van  Dyke  ,  p.  91  and  pp  134—136,  and 
Schlichting  pp  144—  147).  While  higher  order  boundary  layer  theory  does  not  lend  itself 
readily  to  numerical  solution,  it  does  give  a  guide  to  understanding  the  interaction  between  the 
shear  and  core  flows.  Also,  second  order  boundary  layer  theory  contains  all  of  the  terms 
which  are  of  concern  here. 


As  previously  noted,  for  numerical  applications  the  procedure  most  commonly  used  to 
include  the  effects  of  curvature  is  to  include  all  of  the  second  order  terms  in  the  equation  set 
which  is  solved.  Examples  of  this  procedure  are  given  in  References  11,  12,  and  15.  Also 
included  in  Reference  12,  Cebeci  et.al.,  are  curvature  effects  on  the  eddy— viscosity  model. 

The  correction  to  the  outer  inviscid  MOC  solution  is  done  by  displacing  the  wall  by  the 
displacement  thickness,  5*,  and  rerunning  the  inviscid  core  solution.  Examination  of  higher 
order  boundary  layer  theory  shows  that  this  procedure  is  correct  through  second  order.  Hence, 
the  deficiency  with  the  current  boundary  layer  matching  conditions  just  concerns  modifications 
to  the  boundary  layer  codes. 

Wall  Roughness  Effects 

Most  nozzle  boundary  layer  computer  programs  calculate  wall  roughness  effects  in  an 
uncoupled  manner.  For  example,  surface  roughness  height  is  used  to  compute  an  augmented 
heat  transfer  rate.  Such  an  approach  is  completely  unsatisfactory  if  a  purpose  of  the 
calculation  is  to  estimate  nozzle  thrust  loss,  because  the  thrust  loss  calculation  requires  a 
rigorous  energy  balance.  The  wall  roughness  height  found  on  liquid  rocket  wall  materials  are 
small  enough  that  a  coupled  calculation  can  be  performed.  A  coupled  method  for  analyzing 
the  wall  surface  roughness  is  described  below. 


The  effects  of  a  rough  wall  on  the  boundary  layer  can  be  simulated  by  the  method  of 
Cebeci  that  is  described  in  Reference  17.  In  the  method  developed  in  Reference  17,,  the  inner 
region  of  the  Cebeci-Smith  eddy  viscosity  formulation  is  modified  to  provide  for  surface 
roughness.  This  is  done  by  modifying  the  mixing  length 
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I  =  0.4  y  [1-exp  (-y/A)] 


as  described  below.  First,  it  is  recognized  that  the  velocity  profiles  for  smooth  and  rough  walls 
are  similar,  provided  that  the  coordinates  are  displaced,  i.e.,: 

t  -  0.4  (y+Ay)  {  l-exp[-(y+Ay)/A]  } 

The  displacement,  Ay,  is  expressed  as  a  function  of  an  equivalent  sand-grain  roughness 
parameter,  k  ,  i.e.,: 

o 

Ay  =  0.9  (v/U^  [\k+s  -  k+  exp  (-k+/6)] 

where 

k+  =  ks  iyv. 

This  expression  is  valid  for 
4.535  <  k+  <  2000 

where  the  lower  limit  corresponds  to  the  upper  bound  for  a  hydraulically  smooth  surface. 
Typical  values  for  equivalent  sand  roughness,  k  ,  are  given  in  Table  1,  which  has  been  taken 

from  Reference  4. 
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Table  1:  Equivalent  Sand  Roughness  for  Several  Types  of  Surfaces 


Type  of  surface 


ks,cm 


Aerodynamically  smooth  surface 

0 

■7 

Polished  metal  or  wood 

0.05-0.2 

x  10 

Natural  sheet  metal 

0.4 

x  10"3 

Smooth  matte  paint,  carefully  applied 

0.6 

x  10”3 

Standard  camouflage  paint,  average  application 

1 

x  10 

Camouflage  paint,  mass-production  spray 

3 

x  10" 

Dip-galvanized  metal  surface 

15 

x  10~3 

Natural  surface  of  cast  iron 

25 

x  10‘3 

The  metallic  walls,  such  as  Columbium,  of  liquid  rocket  engines  fall  into  the  range  of 

k  given  in  the  above  table.  For  OTV  engine  wall  boundary  layers  the  k*  parameter  should  be 
b  s 

valid,  i.e.,  <2000. 

In  addition  to  modifying  the  eddy  viscosity  model  as  described  above,  it  is  necessary  to 

discuss  the  initial  boundary  layer  profile,  and  the  wall  boundary  conditions.  With  respect  to 
17 

the  former,  Cebeci  used  a  modified  version  of  Cole's  velocity  profile.  This  method  is 

somewhat  complicated  and  is  perhaps  unnecessary.  A  transition  from  a  smooth  wall  profile  by 

the  gradual  introduction  of  roughness  should  allow  the  necessary  profiles  to  be  generated 

automatically.  With  respect  to  the  wall  boundary  condition,  Cebeci  uses  a  displaced  wall 
17 

boundary  condition  .  However,  this  is  not  necessary  and  no  modifications  to  the  wall 
boundary  conditions  were  made. 
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Section  2.2  Implementation 


The  recommendations  of  Section  2.1  were  incorporated  into  the  TDK  code*.  In 
addition,  it  was  discovered  that  the  standard  JANNAF  method  of  computing  the  boundary 
layer  loss  was  inadequate  for  thick  boundary  layers. 

The  following  lists  the  major  changes  to  the  TDK  code.  The  section  numbers  refer  to 
the  sections  in  AFAL  TR— 87-OS  1,  Ref  19  which  is  the  report  which  details  the  TDK  version 
produced  under  Phase  1. 

1)  The  boundary  layer  module,  BLM,  has  been  modified  to  include  the  effects  of 
longitudinal  curvature  on  both  the  mean  flow  and  turbulent  shear  stress  in  the  boundary 
layer.  The  thrust  loss  expression  has  been  completely  re-derived  and  three  options  for 
computing  the  thrust  loss  are  included  in  the  code. 

2)  Section  2.6,  the  BLM  analysis  section,  has  been  rewritten  to  include  the  above 
modifications  and  to  expand  and  clarify  the  numerical  procedure. 

3)  Section  4,  the  program  structure  documentation  has  been  redone  to  reflect  the  code 
changes. 

4)  Section  5.6,  the  BLM  Subroutine  descriptions,  has  been  extended  to  include  the  new 
modifications.  Existing  documentation  has  been  expanded  so  that  the  documentation  is 
more  useful. 

5)  Section  6.8,  the  BLM  Input  description,  has  been  revised  to  reflect  the  changes  in  the 
code. 

6)  Section  7,  the  sample  cases,  has  been  redone  to  reflect  the  output  for  the  current 
version  of  the  code. 

The  TDK  Computer  Program  is  designed  for  engineering  use.  The  FORTRAN  IY 
programming  language  has  been  used  in  an  attempt  to  make  the  computer  program  as  machine 
independent  as  possible.  The  complete  engineering  and  programming  description  of  the  TDK 
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Computer  Program  is  contained  in  Ref.  19. 


Section  2  of  Ref.  19  contains  a  description  of  the  methods  of  analysis  used  in  the 
computer  program. 

Section  3  contains  a  description  of  the  numerical  methods  used  to  integrate  the  fluid 
dynamic  and  chemical  relaxation  equations  in  the  computer  program. 

Section  4  contains  a  description  of  the  program  structure. 

Section  5  contains  a  detailed  engineering  and  programming  description  of  the  program 
subroutines. 

Section  2.3  Phase  I.  Results  and  Conclusions 

Perhaps  the  most  important  result  from  the  Phase  I  work  was  the  discovery  that  the 
standard  JANNAF  performance  methodology  was  inadequate  for  thick  boundary  layers. 

20  21 

The  conventional  JANNAF  method  ’  for  evaluation  of  the  viscous  thrust  loss  in 

rocket  engine  nozzle  employs  both  the  boundary  layer  momentum  and  displacement  thickness. 

22 

This  method  can  be  .jrived  at  by  two  different  approaches  as  shown  by  Alber  .  The  basis  for 
the  approach  is  to  compare  the  thrust  of  an  inviscid  nozzle  to  a  viscous  nozzle  with  the  same 
mass  flow  rate.  In  the  equation  given  below,  the  first  term  represents  the  pressure  forces 
acting  on  an  inviscid  nozzle  and  the  second  and  third  terms  represent  the  total  stress  forces 
acting  on  the  viscous  nozzle. 

AT  -  J  (p*  •  1  <p»  -  pJ2»(rp+S»co4)  dx  (I) 

s 

+  J  Tw  27t(rp+5*cos(j))cos(f)  dx 
o 


whpre  \<z  flip  Hictanrp  frnm  thp  pyic  tn  natAntipl  tupll  onrl  rv*  ic  nr^°.c.!irp  rit 

1 1— V  -  p  IJ  M*«  «»w***S*vw  MtV  |/VkVIM»Uk  TT  U1 1  \%  >  »V4  i>  Lt J  IliV  U1  vJUUI  V  V«1  ip) 


shown  in  Figure  3. 


oc  i  o 

UJ  lO 
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Figure  3.  Nozzle  Wall  Schematic 
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When  the  boundary  layer  characteristics  lengths,  8,  6,  and  $*,  are  assumed  to  be  small 

compared  to  the  characteristic  nozzle  lengths,  r  and  Rc,  then  equation  (1)  reduces  to  the 

21 

standard  JANNAF  relation 


AT  =  2n  costj)  Peu2ro0  -  r0<5*(pc  -  po 


(2) 


If  the  boundary  layer  thickness  assumption  that  .8/ Rc  «  1  is  relaxed,  then  .the  modified 
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JANNAF  relation  can  be  derived  ,  i.e., 


A  i-p  _  — 

=  Pculro0cos(j>  -  r05*(pc  -  pjcostj) 


(3) 


using  the  alternate  forms  of  the  momentum  and  displacement  thicknesses  for  axisymmetric 
flows  shown  below 


(momentum  thickness) 

5  „fl_  -B-  [l  -2-1  dy 
’  ro  Pe.u„  l  u„J  ' 

(4a) 

(displacement  thickness) 

P-hl'-TpST}* 

or0  l  PcucJ 

(4b) 

Table  2  (from  Reference  19)  shows  that  the  differences  in  boundary  layer  loss  as 
computed  from  Equation  1,  2,  and  3  (columns  E,  A,  .and  B  respectively)  are  significant  for 
nozzles  with  thick  boundary  layers. 
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Table  2?  Boundary  Lay'  r  Thru.t  D.ficlt  f or  V.rl.ty  of  Entlnoa  in  tha  Abt.nc. 

and  Fra.anc.  o t  tha  Lontltudlnal  Curvatura.  Cold  Hall  Cut 


THRUST  DEFICIT,  lbf.  aac./lba 


NOZZLE  HITHOUT  LONGITUDINAL  CURVATURE  WITH  LONGITUDINAL  CURVATURE 


EQ.  A 

EQ.  » 

EQ.  C 

EQ.  D 

EQ.  E 

EQ.  A 

EQ.  1 

EQ.  C 

EQ.  D 

EQ.  t 

SSHE 

6.1843 

6.7730 

6.7630 

6.2066 

6.2233 

6.3720 

6.3164 

6.4736 

6.0367 

6.5036 

ASE 

14.1983 

13.6977 

13.4372 

12.3331 

11.9394 

14.2331 

13.7440 

13.3000 

12.3762 

11.7659 

W.-10 

11.6349 

11.0671 

10.8094 

10.4630 

10.4930 

11.7732 

11.0714 

10.6391 

10.4317 

10.4964 

XDELTA 

1.6439 

1.8160 

1.7902 

1.7220 

1.6132 

1.6439 

1.6160 

1.7870 

1.7213 

1.6211 

8C4313 

0.6618 

0.8334 

0.6461 

0.6142 

0.6373 

0.6686 

0.6534 

0.6221 

0.6141 

0.6376 

BC1010 

1.6376 

1.7941 

1.7666 

1.6000 

1.6790 

1.6362 

1.7942 

1.6121 

1.6001 

1.6768 

IUS 

1.9864 

1.9323 

1.9223 

1.8381 

2.0943 

1.9662 

1.9521 

1.6324 

1.6360 

2.0626 

XU-134 

27.3262 

23.6330 

23.3432 

24.6376 

23.2334 

27.3325 

23.6339 

23.2331 

24.6322 

23.1904 

Eq.  A.  AT  =  2n  cos<p  [pg  rQ  0  -  rQ  5*(pe  -  pj] 

Eq.  B.  AT  =  2n  [pg  r0  0  cos <p  -  rQ  5*(pe  -  pj  cos< p) 

Eq.  C.  AT  =  2n  [pe  ujJ  rp  0  cos <p  -  rp  2*(pe  -  pj  cos0] 

Eq.  D.  AT  =  2 n  [p&  u\  rp  0  cosp  -  rp  <5*(pe  -  p  J  cos <p 

‘  ^  pe  ue  ro  ®  cos^  dK  ^  dx] 

O  0 

Eq.  E.  AT  =  1  (p i  *  -  pJ2n  r -  f  (pw  -  PM)2 n  ^  (rS*  cos^dx 

0  F  o 

s 

+  }  tw  2n  (r  +  S*  cos <p)  cos <p  dx 
o  " 

in  all  the  above  equations  e(x)  =  0 
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As  can  be  seen  from  the  above  table,  the  method  by  which  the  boundary  layer  loss  is 
computed  is  more  important  than  the  effect  of  longitudinal  curvatures.  Thus,  this  effect  can  be 
ignored  for  nozzles  of  current  interest. 

Appendix  A  contains  3  published  papers  which  describe  different  facets  of  the  work 
performed  under  Phase  I. 
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3.0  Phase  2  Work  Effort  -  Thick  Boundary  Laver  Assessment 


The  purpose  of  the  Phase  2  '  ork  Effort  was  to  establish,  a)  when  the  thin  shear  layer 
(TSL)  or  boundary  layer  approximations  would  break  down,  b)  the  magnitude  of  the 
breakdown  in  terms  of  performance,  and  c)  to  provide  £  'idance  to  the  new  model  development 
done  in  the  Phase  III  work  effort. 

In  order  to  determine  when  the  TSL  approximation  starts  to  break  down,  a  comparison 
was  made  between  a  traditional  Euler  Solver/Boundary  Layer  Code  (TDK/BLM)23,  a  full 
Navier-Stokes  solver  (VNAP2)24,  and  a  parabolized  Navier  Stokes  solver  (AXI2DS)25.  The 
set  of  engines  examined  by  these  codes  is  given  in  Table  3.  Also  shown  in  this  table  is  the 
relative  boundary  layer  thickness  for  each  of  the  nozzles. 


Table  3:  Boundary  Layer  Thickness  at  the  Exit  Plane  for 
Four  Nozzles.  Cooled  Wall  T  =  1000°R. 

w 


Nozzle 

Area 

Ratio 

Transverse 

Radius  at  the 

Exit  Plane 

r  exit,  inches 

Boundary 

Layer 

Thickness 

8 inches 

5.995/r 

SSME 

77.5:1 

45.12 

1.80 

4% 

204:1 

36.80 

4.40 

12% 

ASE 

400:1 

25.10 

3.25 

13% 

XLR-134 

767:1 

10.97 

2.08 

19% 

Hughes  300:1 

300:1 

1.61 

0.45 

28% 
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The  results  from  this  comparison  are  reported  in  AL-TR-90-041,  the  Phase  >TI  final 
report26,  and  in  Ref.  27,  (which  also  appears  in  Appendix  B).  These  results  were  far  from 
conclusive.  Perhaps  the  greatest  problem  was  the  disparity  in  accuracy  between  the  codes  at 
the  resolutions  we  could  afford  to  use  them.  Since  the  codes  were  the  best  we  could  obtain 
after  a  very  significant  amount  of  effort,  we  concluded  that  both  the  NS  and  PNS  solver  were 
not  up  to  the  tasks  presented  them. 

However,  general  trends  were  observed.  These  trends  were 

1)  When  the  boundary  layer  thickness  became  larger  than  10%  of  the  nozzle,  .the 
departure  between  the  NS,  PNS  solutions  and  the  JEuler/BL  solutions 
accelerates. 

2)  When  compared  to  data,  the  Euler/BL  solutions  are  remarkably  .accurate  even 
for  very  thick  boundary  layers. 

3)  In  terms  of  non-performance  related  parameters,  'boundary  layers  thicker  than 
Sir  >  .1~.,15  are  suspect. 

4)  The  subsonic  portion  of  the  'boundary  layer  is  Infinitesimal  '('^subsonic  <  5/1000) 
-in  the  supersonic  portion  of  large  propulsive  .nozzles. 

With  respect  to  objective  c),  the  conclusions  reached  were  that  a  PNS  solver  was 
desired  for  the  Phase  III  work.  The  reasons  for  this  decision  were 

1)  That  a  PNS  solver  was,  in  terms  of  computer  resources,  the  only  economically 
viable  way  that  both  the  viscous  layer  and  invicid  core  flow  could  .be  .accurately 
resolved. 

2)  That  coupling  a  PNS  solver  to  a  method  of  characteristics  solution  was  not 
desirable  due  to  the  differences  in  numerical  accuracy  .between  ‘these  two 
schemes.  This  rejected  approach  was  what  had  been  originally  proposed. 

As  previously  mentioned,  Ref.  27,  the  technical  paper  reporting  'this  work  is  in 
Appendix  B. 
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4.0 


Phase  3.0  Effort  -  New  Model  Development 


The  purpose  of  the  Phase  3  effort  was  to  develop  a  new  computer  model  which  would 
extend  the  state-of-the-art  in  performance  prediction  for  high  area  ratio  liquid  propellant  rocket 
engine  nozzles.  To  this  end,  a  computer  code  called  VIPER  (Viscous  Interaction  Performance 
Evaluation  Routine)  was  developed.  This  work  is  completely  documented  in  AL-TR-90-042. 
Two  technical  papers  describing  this  work  are  given  in  Appendix  C. 

Following  the  recommendations  of  the  Phase  2  Work.  The  VIPER  was  based  on  the 

Parabolized  Navier  Stokes  equations.  The  PNS  equations  are  a  subset  of  the  Navier-Stokes 

28  29 

(NS)  equations  which  are  valid  for  supersonic  flows  ’  .  The  PNS  equations  neglect  the 

streamwise  diffusion  term  which,  along  with  special  treatment  of  the  subsonic  region  of  the 

boundary  layer,  removes  the  spatial  ellipticity  from  the  steady  form  of  the  equations  and 

permits  a  solution  using  a  streamwise  marching  computational  technique.  Although  the  PNS 

30  31  32 

models  were  developed  in  the  early  1960'so  ’  ’  ,  they  were  not  widely  used  until  the  1970’s 
33-38 

and  1980's 

30  31 

In  the  early  1960's  the  first  PNS  models  were  developed  by  Ferri  and  Morretti  ,  and 
32 

Edelman  .  They  introduced  their  explicit-based  numerical  techniques  to  treat  coupled 
wave-mixing-chemical  kinetic  processes. 

The  PNS  equations  are  integrated  throughout  the  viscous  and  inviscid  regions  of  the 

flow.  This  procedure  eliminates  the  need  to  specify  the  edge  conditions  in  matching  boundary 

layer  and  inviscid  solutions,  i.e.,  the  conventional  inviscid-viscous  interaction.  Popular 
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algorithms  for  solving  PNS  equations  are  those  by  Briley  and  McDonald  and  Beam  and 
40 

Warming  . 

The  VIPER  code  is  a  PNS  solver  for  supersonic  nozzle  flow  with  finite  rate  chemistry. 
The  Beam-Warming  scheme  has  been  employed  to  solve  the  governing  equations.  Richardson 
Extrapolation  for  automatic  step  size  control  in  the  marching  direction  has  also  been 
incorporated  into  the  code.  This  code  can  treat  both  laminar  and  turbulent  flows. 


The  VIPER  code  is  fully  documented  in  AL-TR-90-042,  the  Phase1 3  Report,.  Ref.  41. 
The  code’  was  verified  by  running  7  test  cases.  The  characteristics  of  these  nozzles  are  given 
in  the  table  below. 


Table  4:  Nozzle  Characteristics  for  VIPER  Verification'  Test  Cases- 


Engine’  Name 

Chamber 

Pressure 

Throat 

Radius(in) 

Expansion 

Ratio 

Throat 
Reynolds 
Number  Re] 

NASA/LeRC  Hi-E 

360 

0.5 

1025. 

1.73xl05 

XLR-134 

510 

0.396 

767.9 

1.80x10s 

STS/RCS 

150 

1.021 

28.46 

1.75x10s 

SSME 

3285 

5.1527 

77.5 

1.18xl07 

RL  10 

394.3 

2.57 

205.03 

- 

ASE 

2287 

1.254 

400.7 

2.20x1 06 

HAC-5  lbf 

25 

0.0935 

296.6 

l.lOxlO4 

Cases  1,2,5,  and  6  were  selected  because  they  represent  high  area  ratio  nozzles  which 
are  the  design  point  of  the  VIPER.  Cases  3  and  7  are  small  thrusters  which  have  potentially 
large  boundary  layers.  Lastly,  case  4,  which  is  the  SSME,  was  chosen  to  test  the  high 
Reynold's  number  capability  of  the  code.  Results  for  these  engines  are  tabulated  in  Table  5. 
The  agreement  to  measured  data  (when  available)  is  excellent  when  it  is  realized  that  the 
effect  of  combustion  or  energy  release  efficiency  is  not  included  in  these  calculations.  For 
example,  for  the  NASA/LeRE  hi-expansion  ratio  engine,  the  measured  value  of  C*  efficiency 
was  95.5%  (reading  112  of  Reference  42).  Applying  this  efficiency  to  the  VIPER  prediction 

42 

gives  a  value  of  466.49  lbf-sec/lbm.  Compared  with  the  measured  value  of  468.9  the 
difference  is  -0.51%. 
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Table  5:  Engine  Performance 


Engine  Name 

VIPER 

I  Prediction 
sp 

lbf  -  sec/lbm 

TDK/BLM 

I  Prediction 
sp 

lbf  -  sec/lbm 

Measured 

!sp 

lbf  -  sec/lbm 

NASA/LeRC  Hi-E 

488.47 

(466.49)* 

480.3 11 
(458.70)* 

468.942 

XLR-134 

462.29 

468.681 

- 

STS/RCS 

311.49 

- 

- 

SSME 

455.15 

457.71 

452.61 

RL  10 

462.29 

463.03 1 

458.742 

ASE 

470.20 

473.581 

477.943 

HAC-5  lbf 

218.05 

216.65 

214.5244 

*  corrected  for  95.5%  measured  c*  efficiency 

Comparisons  with  the  JANNAF  standard  performance  prediction  code,  TDK1  are  also 
shown  in  Table  5.  The  comparison  between  the  two  codes  are  best  for  the  high  Reynolds 
number  cases.  The  lack  of  a  consistent  trend  between  VIPER  and  TDK/BLM  is  attributed  to 
the  fact  that  neither  code  includes  combustion  efficiency  correlations. 

Besides  the  good  agreement  between  VIPER  and  TDK/BLM-data,  another  good 
accuracy  check  is  how  well  a  code  conserves  mass,  momentum,  and  energy.  For  finite 
difference  codes  such  as  VIPER,  conservation  of  these  quantities  is  a  very  good  measure  of  the 
overall  validity  of  the  computational  solution  (it,  of  course,  does  not  check  the  validity  of  the 
models  used).  Table  6  gives  a  summary  of  the  conservation  checks  computed  by  the  code. 
Mass  and  momentum  are  conserved  veiy  well  in  all  cases.  Energy  conservation  is  a  much 
more  difficult  quantity  to  achieve.  Partly,  this  is  because  it  requires  good  mass  and 
momentum  conservation,  and  partly  because  the  reference  base  for  energy  is  arbitrary.  As  a 
case  in  point,  the  energy  flow  m(h+V2/2)  of  the  E1000  nozzle  is  only  -0.011  BTU/sec  across 
the  start  line  because  the  chemical  enthalpy  (negative)  balances  out  the  kinetic  energy.  Hence 
a  small  change  in  total  energy  can  result  in  a  very  large  percentage  difference  in  energy 
conservation. 
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Table  6.  Conservation  Checks 


Nozzle 

%  Mass 

Flow  Drift 

%  Momentum 
Flow  Drift 

%  Energy 

Flow  Drift 

NASA/LeRC  Hi-E 

-1.28 

-1.62 

very  high 
near  zero-base. 

XLR-134 

0.02 

1.01 

-6.3 

STS/RCS 

-  0.014 

0.108 

2:9 

SSME 

0.0126 

0.72 

-0.43 

RL  10 

-  0.085 

0.57 

-1.67 

ASE 

-1.18 

-0.31 

-5.3 

HAC-5  lbf 

-0.88 

-1.34 

0.76 

As  is  typical  of  codes  using,  the  Beam-Warming:  difference:  technique,  damping-  of, 
oscillations  caused  by  compression  or  shock  waves  is.  required.  The?  VIPER  code.  uses. explicit 
4th  order  damping,  to  get  rid.  of  these  high  frequency  oscillations.  The.  damping  term  can  be 
written  as 


-  ^  [Q*+2  -  4  Qi+1  +  6  -4  Q^'1  +  O’"2] 

where  is  an  input  number  to  the  code.  A  value  of  =  50  is  recommended.  To  illustrate 

the  effect  the  damping  on  the  flowfield  and  computed  I  ,  a,  study  was  conducted  for  all  7  of 

sp 

the  validation  engines.  Values  of  of  0  and  50  were  used  for  the  comparison.  The  effect  of 

damping  was  dramatic  on  the  flowfield  in  terms  of  shock,  location,  (seeFigurc.4),  however,  the 
difference  in  I  was  less  than  .07%  for  all  cases,  which  indicates  that  the  damping  has  very 

little  effect  on  the  integrated  values  of  wall  pressure  and-.shear  stress.  These  results  are  shown 
in  Table  7.  Figure  5  shows  a  comparison  of  computed  wall  shear  stress,  for*  the.  SSME.  nozzle- 
with  and  without  damping.  Again  the  effects  are  dramatic.  However,  the  difference  in 
integrated  wall  shear  stress  was  only  17  lbf  out  of  a  total  of  7824.  ibf. 
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(  PSF  ) 


Table  7  Effect  of  4th  Order  Damping 
I  (lbf  -  sec/lbm) 


4th  Order 
Damping 

NASA/Lewis  Hi-E 

488.47 

XLR-134 

462.29 

STS/RCS 

311.49 

SSME 

455.15 

RL  10 

462.22 

ASE 

470.20 

HAC-5  lbf 

218.05 

AIM 

304.57 

No 

4th  Order 
Damping 

AIsp 

% 

Difference 

488.39 

0.09 

0.018 

462.00 

0.29 

0.063 

311.44 

0.05 

0.016 

455.18 

0.03 

0.007 

461.99 

0.21 

0.045 

470.19 

0.01 

0.002 

218.04 

0.01 

0.005 

304.64 

0.07 

0.023 

As  with  any  code  completing  its  first  state  of  development,  there  are  deficiencies  and 
unexplained  characteristics.  The  major  problem  areas  are  associated  with  the  grid  generation, 
step  size  control,  the  k-£  turbulence  model,  and  the  sublayer  model. 

The  deficiency  with  the  mesh  generation  and  step  size  affects  the  robustness  of  the 
code.  Values  of  the  stretching  coefficient,  a,  must  be  selected  a  priori  so  that  there  is 
sufficient  mesh  to  resolve  the  boundary  layer  for  the  entire  flow  field.  The  selection  of  the 
marching  step  size  is  also  critical.  It  can  neither  be  too  big  nor  too  small.  The  step  size  can  be 
dynamically  controlled  in  the  current  version  of  VIPER.  However,  there  is  a  distinct  coupling 
between  the  tj  mesh  spacing  and  the  A£  step  size.  The  recommendation  here  is  to  dynamically 
modify  the  mesh  spacing  so  as  to  automatically  resolve  the  wall  shear  layer.  The  step  size 
control  logic  will  then  have  to  be  revised  so  as  to  properly  interact  with  the  dynamic  mesh 
spacing. 

In  liquid  rocket  engines,  the  majority  of  the  vorticity  generated  in  the  flow  comes  from 
the  wall  shear  layer.  Other  sources  of  vorticiiy  are  the  finite  rate  chemistry,  O/F  striations,  and 
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other  phenomena  which  result  in  total  enthalpy  variation.  Because  the  wall  generates  the 
majority  of  the  vorticity  in  the  flow,  the  Cebeci-Smith  eddy  viscosity  model  used  in  VIPER  is 
quite  adequate.  However,  there  are  situations  when  it  would  be  desirable  to  use  a  k-e 
turbulence  model.  The  current  k-e  model  in  VIPER  is  not  robust  and  has  an  inadequate  wall 
treatment.  Hence,  we  recommend  that  the  current  model  be  modified  to  generate  better  initial 
profiles  and  that  the  wall  sublayer  treatment  be  modified  to  incorporate  the  latest  models. 

The  VIPER  code  currently  uses  the  Vigneron  sublayer  model  to  suppress  "departure 
solutions"  in  the  subsonic  portion  of  the  boundary  layer.  As  part  of  that  model,  there  is  a 
safety  factor,  a _,  which  controls  where  the  sublayer  model  is  turned  off.  Numerical 

experiments  have  shown  that  the  results  are  quite  sensitive  to  this  parameter  in  contradiction  to 
most  of  the  literature.  We  recommend  that  the  sublayer  model  be  further  investigated  and 
either  replaced  or  modified  to  control  the  sensitivity. 

Further  improvements  to  VIPER  should  include  the  addition  of  a  discrete  particulate 
phase.  Such  an  addition  would  allow  VIPER  to  analyze  both  solid  and  liquid  propellant 
systems.  The  following  is  a  list  of  potential  particulate  micro-models  which  could  be  included 
in  the  code. 


Micro  Model 

Importance  for 

Liquid  Solid 

Propellant  Propellant 

Vaporizing  Droplets 

V 

Phase  Change 

y] 

Droplet  Super  Cooling 

V 

Breakup 

V 

y/ 

Aglomeration 

V 

V 

Very  Small  Particles  (dp  <  lji) 

yj 

Condensation/Mass  Transfer 
between  Phases 

V 

30 


We  would  recommend  that  all  of  the  above,  except  mass  transfer  between  phases,  be 
incorporated  into  VIPER.  The  basic  particulate  model  would  be  a  fully  coupled  model 
between  the  gas  and  discrete  phase. 

Once  particle  capability  has  been  added  to  the  VIPER,  we  recommend  that  a  simple 
combustion  chamber/motor  cavity  model  be  added.  Such  a  model  would  allow  the  user  to  run 
parameteric  studies  over  a  rational  parameter  space  in  order  to  assess  the  impact  of  each  item 
of  interest. 

In  conclusion,  a  computer  program  called  VIPER  has  been  developed  which  accurately 
predicts  the  flowfield  in  liquid  propellant  rocket  engine  nozzles.  The  VIPER  code  has  been 
developed  as  a  natural  extension  to  the  JANNAF  standard  nozzle  performance  code,  TDK.  By 
incorporating  elements  of  the  TDK  code  in  VIPER,  the  user  community  knowledge  base  has 
been  preserved.  That  is,  users  familiar  with  TDK  can  generate  successful  VIPER  runs  in  a 
matter  of  a  few  hours. 
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5.0 


Ehase  4  Work  Effort  -  Experimental,  Test  Plan 


The  purpose  of  the  Phase  4  effort  was  to  devise  an  experimental  test  plan  which  would 
validate  the  models  developed  in  the  Phase  1  and  3  work  efforts. 

One  of  the  major  problems  with  predicting  the  boundary  layer  loss  in  rocket  exhaust 
nozzles  is  that  there  are  no  direct  measurements  t  '  the  boundary  layer  parameters  in  real 
nozzle  flows  to  use  to  refine  the  analytical  models.  This  observation  is  especially  true  for  the 
turbulence  models  used  in  boundary  layer  calculations.  Current  verification  efforts  are  based 
on  comparison  to  data  in  other  environments,  global  heat  transfer  data  for  real  engines,  or 
backing  out  the  boundary  layer  loss  from  thrust  measurements  by  subtracting  the  other  known 
losses  in  the  system. 

The  other  performance  losses  are  currently  evaluated  by  computer  programs,  such  as 

TDK  ,  BLM  ,  BLIMPJ  ,  etc  using  the  JANNAF  thrust  chamber  evaluation  procedures.  These 

procedures  are  based  upon  a  physical  model  that  accounts  for  the  processes  occurring  in  the 

thrust  chamber,  losses  associated  with  these  processes  and  interactions  among  the 
20  21 

processes  ’  . 

In  the  model,  propellants  enter  the  combustion  chamber  through  the  injector,  are  mixed, 
vaporized  and  combusted.  Deviation  from  complete  and  homogeneous  mixing  vaporization  or 
combustion  to  equilibrium  are  referred  to  as  energy  release  losses  or  injector  losses. 
Theoretically,  these  losses  are  modeled  by  the  Coaxial  Injection  Combustion  Model  (CICM) 

46 

and  the  Standard  Distributed  Energy  Release  Model  (SDER)  codes  .  Experimentally,  the 
energy  release  loss  can  be  estimated  from  the  characteristic  exhaust  velocity  efficiency, 


The  products  of  the  chemical  reactions  are  then  expanded  in  the  nozzle.  The  reactions 
continue  during  the  expansion.  Deviation  from  the  local  chemical  equilibrium  are  referred  to 
as  the  kinetic  losses.  There  are  no  direct  methods  of  measuring  the  kinetic  losses, 


The  losses  due  to  non-unifonn  expansion  of  the  available  momentum  in  the  direction  of 
thrust  are  referred  to  as  two-dimensional  or  divergence  losses.  Theoretically,  these  losses  are 
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(Method  of  Characteristics)  modules  of  the  TDK  code. 
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Graphically,  all  of  the  above  described  losses  are  shown  in  Figure  6. 
performance  is  based  on  ODE  calculations.  _ . 


oot  specific 
impulse 


heat  loss 


PEL  IVERED  THMJST 
CHAH6EA  SPECIFIC 
IMPULSE 


The  ideal 


Figure  6.  Illustration  of  thrust  chamber  losses 
from  ideal  performance 


Systematic  measurements'  to  verify  the'  predicted'  magnitude'  of  the;  aforementioned' 
losses  are  sparse;-  and  a'  high1  degree  of-  uncertainty  is  associated  with’  the  available  data;.  For 
both  the  kinetic  and  boundary  layer  losses,  experimental  validation of the predictive  techniques 
fof  real  engines  are  not-  available.  For  low  area  ratio  nozzles,  the  boundary  layer  loss  has-  been 
considered  small  enough  that-  any  errors  associated  with  the  predictions  were  acceptable.  The 
magnitude  of  the  boundary  layer  loss  as  a  percent  of.  total  per  ormance  increases  with  area 
ratio.  Figure  7  shows-  the'  boundary  layer  losses  as  calculated  by  the  SPP  and  TDK  codes  for  a 
Variety  Of  motors  as  a  function  of  area  ratio.  The  dependency  of  ‘hese  results  on  the  particular 
turbulence'  model  used  in'  the'  calculation  are  on  the  order  of  ±  30-40%  .  That  is,  a  predicted 
10  second  loss  could  actually  be  only  7  seconds  or  could  be  as  high  as  13  seconds.  Since  a  1% 
change  in  Isp  caff  result  iff  approximately  a  4%  change  in  payload,  the  accurate  knowledge  of 
the  boundary  layer  loss  is  very  important.  Hence,  the  boundary  layer  calculations  must  be 
validated  directly  using  experimental  data. 

In  order  to  validate  or  extend  analytical  methods,  the  planned  experimental  effort  must 

include  techniques  to  measure  boundary  layer  parameters  that  will  produce  high  quality  data. 

Qualify  data  on  nozzle  wall  boundary  layers  is  very  limited  at  present.  Perhaps,  the  best  work 
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in  this  area  is  that  of  Back,  et  al,  at  JPL  ’,  which  was  for  cold  flow. 

As  part  of  the  effort  to  prepare  an  experimental  validation  test  plan,  numerous  testing 
facilities  were  contacted  and  visited.  Available  experimental  methods  and  facility  and  test 
requirements  were  reviewed,  and  recommendations  were  made  to  the  Astronautics  Laboratory. 
Results  of  the  study  were  documented  in  the  Phase  4  final  report  together  with  a  test  plan 
intended  to  validate  the  analytical  methods  developed  during  the  boundary  layer  study. 

Certain  ground  rules  were  used  in  preparing  this  test  plan.  The  first  ground  rule  was 
that  we  wanted  a  program  which  would  measure  the  desired  data,  we  did  not  want  to  embark 
upon  a  research  program  to  develop  improved  boundary  layer  diagnostic  techniques. 
Secondly,  we  wanted  the  test  environment  to  simulate  as  close  as  feasible  the  conditions  in  a 
real  propulsive  nozzle.  The  minimum  requirements  for  this  condition  are  shown  in  Table  8, 
below.  Last,  and  most  important,  we  wanted  accurate  data. 
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Table  8.  Miminum  Test  Environment  Criteria 


o  The  flow  in  the  boundary  layer  must  be  turbulent  from  before  the  throat. 

o  The  boundary  layer  should  encompass  between  10%  to  40%  of  the  flow  at  the 
exit  plane. 

o  The  heat  transfer  to  the  wall  should  be  high  enough  to  simulate  conditions  in  a 
real  rocket  nozzle. 

c  The  expansion  process  should  maintain,  as  realistically  as  possible,  the 
conditions  in  a  real  rocket  nozzle.  For  example,  axisymmetric  instead  of  2D 
nozzles,  ratio  of  specific  heats  <  1.4. 

Some  of  the  ground  rules  were  found  to  be  mutually  exclusive.  For  example,  accurate 
data  taking  precluded  using  a  real  engine  as  a  test  bed.  The  environment  in  real  engines  is  too 
severe  for  many  diagnostic  techniques.  We  also  found  that  off  the  shelf  techniques  were  not 
sufficiently  accurate  for  all  of  the  data  measurements  desired. 

Parameters  to  be  Measured 

The  primary  interest  is  computing  the  boundary  layer  parameters  in  propulsive  nozzles 
is  in  determining  the  viscous  thrust  loss  and  wall  heat  transfer  rate.  The  most  fundamental 
property  in  the  thrust  loss  is  the  wall  shear  stress.  The  heat  transfer  to  the  wall  is  determined 
by  the  temperature  (enthalpy  for  reacting  flows)  gradient  at  the  wall.  While  both  the  shear 
stress  and  temperature  gradient  can  be  deduced  from  other  quantities,  a  direct  measurement  of 
these  items  is  definitely  preferred.  In  the  absence  of  direct  measurements,  diagnostics  which 
require  the  least  amount  of  assumptions  are  preferred. 

Average  velocity  and  temperature  throughout  the  boundary  layer  and  Reynolds  stress 
terms,  u',  v',  are  other  parameters  to  be  measured.  Although  the  last  two  are  less  important, 
their  correct  evaluation  yields  valuable  data  for  validation  of  turbulence  models.  Briefly  the 
following  measurements  are  to  be  made: 
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1.  Local  wall  shear; 

2.  Local  heat  flux  at  the  wall; 

3.  Average  velocity  and  temperature  in  the  boundary  layer; 
and 

4.  Reynolds  stress  terms. 

Conclusions 

After  talking  with  various  organizations,  we  concluded  that  two  families  of  techniques 
are  suitable  for  hypersonic  flow  velocity  measurements:  particle  scattering  and  molecular 
scattering.  Particle  scattering  techniques  which  include  real  fringe  laser  Doppler  velocimetry, 
laser  transit  anemometry,  and  Doppler  spectrometry,  are  well  developed  and  velocity 

measurements  in  supersonic  and  hypersonic  flow  regimes  have  been  demonstrated.  However, 

paidcle  scattering  techniques  require  seeding  the  flow,  which  introduces  ambiguity  due  to 
particle  lag  in  accelerated  or  decelerated  flows.  The  question  of  particle  lag  has  been 
addressed  and  it  is  concluded  that  the  measurement  uncertainties  are  small  when  the  particles 
used  are  less  than  0.3  [im  in  shock-free  expansion  flows. 

The  uncertainty  of  the  fluctuating  components  of  the  velocity  within  the  boundary  layer 
will  depend  on  the  density  of  the  flow  and  the  frequency  response  of  the  particles. 

Molecular  techniques  which  encompass  laser  induced  fluorescence  (LIF),  Inverse 
Raman  Scattering  (IRS)  and  Coherent  Anti-Stokes  Raman  Scattering  CARS,  are  being 
developed.  These  techniques  have  the  potential  of  yielding  more  accurate  results  than 
particle-based  techniques.  However,  several  practical  aspects  hinder  their  application  to 
hypersonic  flow  measurements.  LIF  of  alkalis,  which  is  more  accurate  than  LIF  of 

requires  seeding  the  wind  tunnel  with  highly  corrosive  materials  such  as  sodium.  Seed 
molecules  condense  under  the  low  pressure  and  temperature  conditions  characteristic  of 
hypersoi  ic  wind  tunnels.  LIF  of  is  still  under  investigation  and  preliminary  analysis  shows 

that  the  velocity  error  is  about  500  m/sec.  IRS  and  CARS  measure  the  velocity  from  the 
spectrum  of  the  coherent  Raman  scattered  signal,  and  therefore,  pressure  broadening  due  to  a 
shock  wave  can  result  in  large  errors.  Velocity  measurement  using  IRS  is  on-axis,  and  hence 
the  spatial  resolution  is  poor.  IRS  and  CARS,  however,  have  the  advantage  of  yielding  the 
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temperature  and  pressure  simultaneously  with  the  velocity,  giving  a  unified  approach  to  flow 
field  measurement. 

Table  9  is  a  summary  of  velocity  measurement  techniques.  The  table  shows  the 
demonstrated  velocity  range,  accuracy,  SNR  (signal  to  noise  ratio)  and  other  relevant 
parameters  of  techniques. 

Sgg.Q.oim£n.daii.Qii 

The  primary  objective  of  developing  an  experimental  plan  for  boundary  layer 
measurements  is  to  obtain  accurate  and  high  quality  data  to  compare  against  the  present  and 
future  analytical  models.  During  the  course  of  our  investigations,  no  single  or  set  of  diagnostic 
techniques  and  experimental  facilities  were  found  to  be  capable  of  achieving  this  objective. 
As  a  result,  it  was  necessary  to  relax  some  of  the  ground  rules  which  we  had  originally  laid  out 
in  order  to  meet  the  primary  goal  of  obtaining  accurate  data.  As  mentioned  in  the 
introduction,  we  concluded  that  hot  engine  firings  were  not  compatible  with  accurate  mapping 
of  the  boundary  layer.  This  conclusion  was  based  on  considerations  of  available  experimental 
techniques,  test  times,  facility  requirements,  and  cost.  We  also  discovered  that  there  were  no 
diagnostic  techniques  which  were  "off  the  shelf'  ready  for  boundary  layer  measurements  of  the 
required  accuracy  in  simulated  rocket  nozzles.  In  addition,  facility  capabilities  for  different 
organizations  varied  tremendously.  Calspan  has  excellent  short  duration  test  facilities,  while 
AEDC  has  good  large  scale  longer  duration  capability.  However,  the  selected  diagnostic 
techniques  must  be  mated  with  the  correct  facility  to  insure  adequate  data  acquisition. 
Because  of  the  developmental  nature  of  the  diagnostic  techniques  and  the  mating  of  these 
techniques  with  the  facilities,  the  following  recommendations  are  conditional  rather  than 
absolute.  First,  we  will  cover  the  recommended  diagnostic  methods  and  then  discuss 
recommendations  for  test  facilities. 

Diagnostic  Method  Recommendations  for  Velocity  Profile  Measurements 

LDV  must  be  considered  the  primary  candidate  for  taking  velocity  measurements  in  the 
boundary  layer.  This  conclusion  is  based  on  the  fact  that  this  method  is  well  developed  in 
comparison  to  other  laser  techniques.  Alternates  and/or  backup  methods  should  include  LIF 
and  CARS.  We  strongly  recommend  that  optical  techniques  should  be  the  cornerstone  of  these 
measurements. 
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Summary  of  Velocity  Techniques,  Reference  2  5 


We  view  both  LIF  and  CARS  as  potentially  superior  methods  of  velocity  measurement. 
The  trade  off  between  LDV  and  LIF/CARS  is  the  development  time  and  difficulties  associated 
with  particulate  seeding  for  LDV  versus  moving  LIF/CARS  from  the  lab  into  a  test 
environment.  Should  the  latter  be  accomplished,  then  LIF/CARS  would  be  the  recommended 
technique. 

Static  Pressure  Measurements 

Pressure  transducers  on  the  wall  are  the  recommended  diagnostic  technique  for 
measuring  static  pressure  at  the  wall.  If  the  selected  test  conditions  are  benign  enough,  pitot 
measurements  can  also  be  used  as  a  backup  and  to  verify  calibration  of  both  the  wall  static 
pressure  and  velocity  profile. 

Temperature  Profiles 

Laser  based  non-intrusive  techniques  are  recommended.  LIF  seems  to  have  a  slight 
edge  on  CARS  based  techniques,  but  both  methods  should  perform  well  theoretically.  The 
final  choice  should  be  decided  by  the  selected  contractor  based  on  cost  and  risk  factors.  Both 
LIF  and  CARS  are  also  capable  of  determining  the  density  profiles  in  the  flow.  Information 
about  the  chemical  composition  would  thus  supply  enough  information  to  compute  j  ’ 
pressure  profile.  The  redundancy  in  pressure  data  will  supply  a  useful  accuracy  chf 
measurements.  Total  temperature  probes  are  also  readily  available  (see  Appendix  .  . 
should  be  used  as  a  calibration  backup  if  pitot  measurements  are  also  taken. 

Wall  Heat  Transfer  Measurement 

There  are  any  number  of  adequate  methods  to  measure  the  wall  heat  transfer  rate. 
However,  most  of  these  methods  depend  on  the  temperature  range  of  the  flow  (total  and  static) 
and  of  the  wall.  The  final  selection  must  be  based  on  the  test  environment  selected.  However, 
we  do  not  recommend  that  the  heat  transfer  rate  be  deduced  solely  from  backwall  temperature 
measurements.  Such  methods  require  operation  to  steady  state  to  supply  the  desired  accuracy. 


39 


Wall  Shear  Stress 


Direct  measurements  of  tu:  wall  shear  stress  are  definitely  preferred.  If  the  test 
environment  permits,  some  sort  of  the  floating  element  method  should  be  used  (See  Appendix 
C  and  also  the  description  of  the  Kistler  Skin  Friction  Gauge).  Extrapolation  of  the  velocity 
profile  to  the  wall  can  also  be  used  if  floating  elements  are  not  practical.  Under  ho 
circumstances  should  the  Reynolds  Analogy  be  used  as  the  primary  method  of  determining 
skin  friction.  The  use  of  the  law  of  the  wall  also  presupposes  the  form  of  the  results  and  lienee 
its  use  is  discouraged. 

Test  Facilities 

The  ideal  facilities  for  these  tests'  would  be  space  based'  so  that  practical  considerations 
such  as  pumping  requirements,  total’  pressures,,  and  total  temperatures  would1  not  be  of  concern. 
However,  since  we  are  lacking  permanent  space  stations  or  lunar  facilities,  we  must  make  due 
with  earth  based  testing.  It  is  always  preferred’  that-  the  test  environment  closely  simulate  the 
environment  to  which  the  data  is  to  be  applied’.  However,  in  this  case,  we  recommend  against 
real  engine  firings.  Too  many  of  the  recommended*  diagnostic  techniques  are  not  well  enough' 
established  to  allow  for  their  extension*  to- such  a  hostile  environment.  In’ fact,  some  of  the  test 
methods  suggested  would  best  be  done  with  windows  cut  into  the  nozzle  wall.1  Many  of  the 
problems  associated*  with  facility  requirements  disappear  if  the  diagnostic  techniques  can  be 
made  to  work  for  short  duration  test.  Calspan  has  an  excellant  short  duration  test  capability 
which  would  allow  good*  simulation  of  high  area  ratio  nozzles.  However,  it  is  not  clear  if 
either  LIF  or  CARS  methods  can-  be  made  to  work  with  test  durations  in  the  50  millisecond 
range.  The  Accurex  arc  tunnel  suggested  by  one  of  the  experimental  test  plan  contributors 
would  only  be  acceptable  if  the  total  enthalpy  gradient  across  the  nozzle  that  is  produced  by 
the  arc  could  be  eliminated. 

In  order  to  simulate  a  real  nozzle  as  closely  as  possible,  we  make  the*  following 
recommendations. 


1  \ 


The  test  nozzle  should'be  axisyimiictric  and  have  at  least  an  area  ratio  of  300:  i. 


2)  The  flow  should  be  fully  turbulent  before  tun  throat  plane: 
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3) 


The  flow  shouid  be  as  hot  as  is  compatible  with  the  diagnostic  techniques. 
Flows  in  the  range  of  1500°R  are  recommended.  These  types  of  temperatures 
will  allow  for  adequate  heat  transfer  rates,  density  variations,  and  avoid 
condensation  problems. 

4)  The  gas  supply  system  should  supply  nearly  uniform  property  flow. 

In  order  to  assure  that  the  experiments  achieve  their  goal  of  supplying  data  to  validate 
analytical  models,  we  also  recommend  that  an  organization  such  as  Software  and  Engineering 
Associates,  Inc.  be  assigned  to  monitor  the  experimental  work  effort. 
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NOMENCLATURE 


c 

speed  of  light 

dp 

droplet  radius 

g(v) 

line  shape  function 

h 

enthalpy 

I 

intensity 

M 

Mach  number 

n 

normal  direction 

P 

pressure 

Pr 

Prandtl  number 

r 

radius 

★ 

r 

throat  radius 

Re 

Reynold’s  number 

s 

entropy 

T 

temperature 

u 

velocity 

U 

velocity 

V 

velocity 

V 

velocity  (total) 

y 

distance 

Greek 


9  boundary  layer  momentum  thickness 

c,  vorticity 

p  density 

8  boundary  layer  thickness 

A  difference 

t  shear  stress,  characteristic  time 

H  coefficient  of  viscosity 

v  frequency 
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NOMENCLATURE  (Continued) 


Abbreviations 


AL 

ASE 

BC4515 

BC1010 

BLM 

CARS 

CRC 

IRS 

IUS 

JANNAF 

LDA 

LDS 

LDV 

LIF 

lox/gh2 

LTA 

RL-10 

SEA 

SRGS 

SSME 

TDK 

XDELTA 

XLR134 


Astronautics  Laboratory 

Rocketdyne  Advance  Space  Engine,  area  ratio  400: 1 

Experimental  Nozzles  by  JPL,  cone  shaped 

Boundary  Layer  Module  of  TDK  computer  code 

Coherent  Anti-Stokes  Raman  Spectroscopy 

Coherent  Raman  Spectroscopy 

Inverse  Raman  Spectroscopy 

Inertial  Upper  Stage,  Space  motor,  OTV 

Joint  Army  Navy  NASA  Air  Force 

Laser  Doppler  Anemonetry 

Laser  Doppler  Spectrometry 

Laser  Doppler  Velocimetry 

Laser  Induced  Fluorescence 

Liquid  Oxygen  and  Gaseous  Hydrogen 

Laser  Transit  Anemonetry 

Pratt  &  Whitney  Space  Engine,  area  ratio  205:1 

Software  and  Engineering  Associates,  Inc. 

Simulated  Raman  Spectroscopy 
Space  Shuttle  Main  Engine,  area  ratio  76:1 
Two  Dimensional  Kinetic  Computer  Code 
Extended  Delta  (Solid  Propellant  Space  motor) 

OTV  Engine,  Space  Storable,  area  ratio  767:1 
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ABSTRACT 


The  performance  loss  in  high 
expansion  ratio  nozzles  due  to  tne 
turbulent  boundary  layer  is  a  critical 
design  parameter  in  selecting  the  nozzle 
geometry.  This  wor;<  studies  the  effects 
of  transverse  and  longitudinal  curvature 
on  the  boundary  layer  growth  and  their 
Impact  on  performance.  Curvature  effects 
become  important  in  axisymmet r ic  flows 
when  the  curvature  radius  is  of  the  same 
order  as  the  boundary  layer  thickness. 

The  effects  of  longitudinal  curvature 
manifest  themselves  as  a  centrifugal 
force  term  on  the  mean  flow  in  the 
boundary  layer,  and  also  on  the  turbulent 
shear  stress  model.  The  present  Boundary 
Layer  Module  ( BLM )  of  the  TDK  Computer 
Code  has  been  modified  for  longitudinal 
curvature  effects  by  the  addition  of 
second  order  boundary  layer  terms  to  the 
equation  set.  The  normal  momentum 
equation  has  also  been  included  in  the 
formulation.  The  edge  conditions  of  the 
boundary  layer  have  been  modified  to 
provide  the  correct  matching  with  the 
outer  flow  expansion.  A  comparison  of 
results  between  the  modified  version  of 
the  code  and  previous  version  of  the  code 
has  been  made.  The  modified  code  has 
been  tested  for  a  variety  of  nozzles  and 
the  results  are  presented. 
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Given  by  Equation  (15a) 
Given  by  Equation  (15b) 
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pressure 
prandtl  number 
Radius  of  longitudinal 
curvature 
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I.  INTRODUCTION 


The  trend  toward  very  large  area 
ratio  nozzles,  which  result  in 
performance  gains  for  space  applications, 
has  Increased  the  need  for  detailed 
knowledge  of  the  momentum  tosses  due  to 
viscous  effects  in  propulsive  nozzles. 

r  *!:ral!lLlondl  aPProachesO-S)  make  use 
or  the  Pr-andM  thin  shear  layer 

approximation  to  compute  the' losses  cue 
to  viscous  affects,  i.e.,  g/t  <<  , 
which  eliminates  the  normal  momentum 
equation  or  the  "centrifugal  force". 
However,  it  should  be  noted  that  the 
normal  pressure  gradient  vanishes  not 
because  the  shear  layer  13  thin,  but 
because  the  surface  1s  relatively  fiat, 
i.e,,  the  radius  of  curvature,  S  ,  is 
iarge  or  more  precisely  R  /%>>?, 

Van  Dyke^  has  shown  that 

longitudinal  curvature  makes  a 


contribution  that  Is  additive  to  that  of 
transverse  curvature  and  of  the  same* 
relative  order.  The  per  t  urb<i  1 1  on  theory 
applied  by  Van  Dyke  demonstrates  the 
Importance  of  second  order1  terms.  The 
metric  influence  of  curvature  has  been 

given  by  Scnu l t z- Cr unow  and  Brower^1 ^ 
for  laminar  and  Incompressible  flows. 
Their  approach  consists  of  adoption  of  a 
curvelinear  coordinate  system  and  taking 
the  arc  length  along  the  surface  as  the 
coordinate.  This  brings  the 
longitudinal  curvature  directly  into  the 
problem  since  any  differentiation  with 
respect  to  x  carries  a  factor  of  the 
ratio  of  curvature  radii  (wall  position 
to  actual  position). 

Cebeci,  Hirsch,  and  Whitelaw^1^ 
analyzed  the  turbulent  boundary  layer  on 
a  (convex)  longitudinally  curved  surface 
using  the  mean  flow  equations  of 

Shu. 1 1 z-Gr unow  and  Breuer^11^  discussed 
earl  \r .  In  their  treatment  they 
replied  the  laminar  viscosity  with  an 
effective  eddy  viscosity.  The  eddy 
viscosity  was  specified  to  be  the 

standard  Cebecl  &  Smlth^^  formulation 

MM) 

modified  by  Bradshaw's  correction  for 

longitudinal  curvature.  Eghllma 
( l  b) 

et.al.  used  a  similar  approach. 

However,  in  ail  these  approaches  the 
assumption  has  been  made  that  the 
vorticity  vanishes  along  the  edge  of  the 
boundary  layer,  i.e.,  the  outer  flow  is 
i rrolat ional . 

Most  of  the  existing  literature  is 
limited  to  boundary  layers  in  external 
flows,  which  have  received  much  more 
attention  than  internal  flows  with  thick 
boundary  layers.  However,  Whitefleld  and 

Lewis^^  have  reported  results  from 
experiments  on  laminar  boundary  layer 
development  in  nozzles.  In  their 
experiments  the  boundary  layer  thickness 
were  as  high  as  90*  of  the  nozzle  radius. 

The  purpose  of  this  work  is  to  study 
the  effect  of  longitudinal  curvature  on 
internal  flows  and  specifically  nozzles. 

To  do  this,  the  problem  has  been 
formulated  and  the  new  terms  In  the 
resulting  donations  in  similarity  form 
have  been  added  to  tho  Bl.M  Comoutor 

(  7  ) 

Code  which  already  accounted  for 
axisymmetric  flows  with  transverse 
curvature.  The  modified  code  has  been 
tested  for  a  variety  of  rocket  engine 
nozzles  and  the  results  have  been 
presented.  The  code  also  has  been  tested 
for  external  flows  and  a  comparison  with 
existing  literature  has  been  made.  For 
"relatl vely"  thick  boundary  layers  in  the 
nozzles,  all  second  order  terms  are 
accounted  for  and  included  in  the 
Boundary  Layer  Module  ( BLM )  of  the  TDK 
Computer  Code.  a  o 


II.  THEORY 


For  a  compressible  boundary  layer 
flow  in  an  uxlsymmetrie  nozzle,  the 
governing  equations  including  tho 
transverse  and  longitudinal  curvature 


effects  are  given  as 
Continuity: 


(12. 15. 17, 18) 
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The  influence  of  longitudinal 
curvature  on  the  mean  boundary  layer  flow 
is  modeled  by  inclusion  of  the 
centrifugal  force  term  in  the  normal 

momentum  equu t i on ^ *  1 0 ' 1 1  V.  (See 
Equation  2b) . 

In  these  equations  n  is  the  flow 
index  which  is  zero  for  two-dimensional 
flow  and  one  for  axisymmetric  flows,  and 
-k"k(x)  Is  the  longitudinal  curvature. 

Note  that  the  -  sign  in  1tk=(x)y  refers 
to  concave  and  ♦  sign  refers  to  convex 
surfaces.  Figure  1  shows  the  definition 
of  r,rQ,  R  and  R . 

Transverse  curvature  is  defined  as: 

y  Cos  $  r 

t  -  - ±  -  L_  { u.i ) 


and  the  longitudinal  curvature  a3: 


-  1  t  k(x)y 


The  streamwise  momentum  equation, 
(2a),  requires  3p/3x  along  lines 
of  constant  y*  From  Sehul tz-Grunow ( \\ ) : 


The  boundary  conditions  for  equations 
(1-3)  are: 


at  y-0,  u-0,  v-vw(x), 


and  T  •  Ty  (x)  or  $w  (x)  is  known;  (5b) 


at  y  -  yfl 


u  -  V  H  “  V 


(lx)y  •  (ix}n  *  ^ly’x  ^x  (lx}y  <*>> 


Surfaces  of  constant  n  are  all 
parallel  to  the  curved  wall.  At  the 

MedgeH  the  term  <!$>  is  given  by 
o  x  n 

Bernoulli's  equation; 
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The  boundary  condition  for  p(y  ) 

/ n  j/j\  ma x 

was  taken  from  Van  Dyke  '  who  matched 
the  "inner  Limit  of  the  Outer  Expansion" 
(J)th  iteration  to  the  (J-I)th  iteration 
of  the  "inner  flow"  (i.e.,  the  boundary 
layer).  Van  Dyke  showed  that  equation 
(2a)  and  (5c)  include  all  second  order 
terms  caused  by  longitudinal  curvature 
effect . 

In  the  governing  equations  the  terms 
for  normal  stresses  have  been  neglected. 
After  Cebeci-Smith  we  define: 


(r*)  (") 

3n  x  3x  y 

Replacing  the  second  term  on  the 
right  hand  side  in  the  above  equation  by 
using  Equations  (2b)  and  (yb)  gives: 


-  p  u  —  £  -  :£*4>~  u2r2 
dx  pe  e  ax  1-kTxJy  e 
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To  obtain  the  equations  in  similarity 
form  a  stream  function  ^  has  been 
defined  such  that: 

„„,.n  -  ( 


The  second  term  on  the  right  hand 
side  is  of  second  order  and  is  included 
in  the  streamwise  momentum  equation  in 
the  mu  coefficient,  m.  <tnd  m-  are  given 
by  Eq?  (-13d).  1  * 


In  terms  of  new  variables  the 
continuity  and  streamwise  momentum 
equations,  l-2a  can  be  rewritten  as: 


p*r"  ■  “kUTy  * (pv )w  ro  ~  1^  ) 


t  ft 
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where  *  represents  and 
on 


and  ;  and  n  are  the  similarity 
coordinates  given  by: 


A  -  Mc>-w-0(ii,c)(r-)  - ----  C  Re 


(i-l)n  ^  (15 a) 
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du  . 

m2  "  u“  dc~  ’  %  *  He  (p.juc) 


tht*  boundary  layer  for  ASE  nozzle  at  the 
exit  plane  *  A  pressure  difference  of 
about  15%  is  observed  across  the  boundary 
layer  with  a  negative  slope  (concave 
surface),  l  .  e  . ,  the  pressure  at  the  wall 
is  15*  higher  than  the  pressure  at  the 
edge  of  boundary  layer. 

For  nozzles  with  “higher"  curvature, 
the  pressure  gradient  across  the  boundary 
layer  may  vary  significantly.  For  SSME 
nozzle  a  pressure  difference  of  about  35* 
Is  observed  from  Figure  6  across  the 
boundary  layer  with  U/R0>exit  ’  °'016, 

The  boundary  layer  is  considerably 
( 22 ) 

thicker  for  SSME^ 

IV  .  CONCLUDING^ REM ARKS 


The  results  in  this  work  indicate 
that  the  longitudinal  curvature  is  of 
importance  in  high  expansion  area  ratio 
and/or  strongly  curved  nozzle  contours. 
The  primary  effect  is  on  the  velocity 
profile  which  results  in  thicker  boundary 
layer.  The  magnitude  of  6 / R  specifies 
the  Importance  of  longitudinal  curvature. 
Also,  it  is  found  out  that  generally  the 
inclusion  of  longitudinal  curvature 
decreases  the  thrust  deficit  and 
consequently  increases  the  performance 
for  rocket  engine  nozzles.  The  method 
used  to  obtain  the  curvature  seems  to  be 
adequate  and  more  stable  and  more 
accurate  than  methods  which  lead  to 
calculation  of  the  second  derivative. 
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Figure  3:  Curvatures  for  ASK  and  SSMK 
Nozzles. 
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Figure  4:  6/R  for  ASE  and  SSMR  Nozzles 
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Across  the  Boundary  Layer  for  ASE  >07.fcie. 
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ABSTRACT 


The  performance  losses  in  rocket  nozzles  due  to  the  viscous  effects  have  been  studied. 
Special  consideration  was  given  to  nozzles  with  high  expansion  ratios.  The  formulation  of  the 
boundary  layer  equations  used  in  this  study  includes  the  effects  of  transverse  and  longitudinal 
curvatures.  The  effects  of  longitudinal  curvature  manifest  themselves  as  a  centrifugal  force 
term  on  the  mean  flow  in  the  boundary  layer,  and  also  on  the  turbulent  shear  stress  model. 
The  results  indicate  that  although  the  longitudinal  curvature  creates  a  fairly  strong  pressure 
difference  across  the  boundary  layer,  the  effects  on  the  perfomiance  is  minimal.  The  results 
also  show  that  in  high  expansion  ratio  nozzles,  the  viscous  layer  becomes  very  "thick"  and  the 
traditional  boundary  layer  assumptions  cause  significant  error  in  the  viscous  loss  calculations. 
Improvements  to  the  method  of  evaluating  thrust  loss  arc  presented. 
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NOMENCLATURE 


H(rj,£)  Enthalpy 

k(x)  or  k(Q  Transformed  longitudinal  curvature,  positive  for  convex  and  negative  for 

concave  surfaces 

L  Length  scale 

p  Pressure 

pr  Prandtl  number 

R  Radius  of  longitudinal  curvature 

r  Radius  of  transverse  curvature 

Re  Reynolds  number 

t  =  y  cos<j>/r0 

T  Thrust 

u,v  Velocity  components 

x,y  Coordinate  system 


Greek  Letters 

5 

8* 

e 

e(x) 

M- 

v 

v 

p 

e 

X 


Boundary  layer  thickness 
Displacement  thickness 
Eddy  viscosity  term 
Given  by  eg. 

Wall  angle 
Viscosity 

Kinematic  viscosity 
Transformed  coordinates 
Stream  function 
Density 

Momentum  thickness 
Shear  stress 
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Superscripts 


* 

d/dr| 

It 

avan2 

* 

Potential  Flow  (at  distance  8*  from  the  wall) 

— 

Modified  Value 

Subscripts 

w 

At  the  wall 

0 

To  the  wall 

u 

Refers  to  velocity 

H 

Refers  to  enthalpy 

c 

At  the  edge  of  boundary  layer 

p 

At  the  edge  of  potential  flow 

Abbreviations 

ASE 

BC4515) 

BClOlOj 

Rocketdyne  Advance  Space  Engine,  area  ratio  400: 1 

Experimental  Nozzles  by  JPL,  cone  shaped 

BLM 

Boundary  Layer  Module  of  TDK  computer  code 

IUS 

Inertial  Upper  Stage,  Space  motor,  OTV 

JANNAF 

Joint  Army  Navy  NASA  Air  Force 

Lox/GH2 

Liquid  Oxygen  and  Gaseous  Hydrogen 

RL-10 

Pratt  &  Whitney  Space  Engine,  area  ratio  205:1 

SSME 

Space  Shuttle  Main  Engine,  area  ratio  76:1 

TDK 

Two  Dimensional  Kinetic  Computer  Code 

XDELTA 

Extended  Delta  (Solid  Propellant  Space  motor) 

XLR134 

OTV  Engine,  Space  storable,  area  ratio  767:1 
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INTRODUCTION 


The  trend  toward  very  large  area  ratio  nozzles,  which  result  in  performance  gains  for 
space  applications,  has  increased  the  need  for  detailed  knowledge  of  the  momentum  losses  due 
to  nozzle  viscous  effects  (i.e.,  boundary  layer)  in  propulsion  systems.  These  losses  degrade 
overall  system  performance  by  increasing  system  weight,  decreasing  useful  payload  weight, 
and/or  decreasing  effective  system  range.  The  traditional  approaches1**  use  the  Prandtl  thin 
shear  layer  approximation  i.e.,  8/L  «  1,  to  compute  the  losses  due  to  viscous  effects.  These 
assumptions  eliminate  the  normal  momentum  equation  or  the  "centrifugal  force  balance".  The 
normal  pressure  gradient  vanishes  not  because  the  shear  layer  is  thin,  but  because  the  surface 
is  relatively  flat,  i.e.,  the  radius  of  curvature  is  large. 

Van  Dyke9-10  has  shown  that  longitudinal  curvature  makes  a  contribution  that  is 
additive  to  that  of  transverse  curvature  and  is  of  the  same  relative  order  of  magnitude.  The 
perturbation  theory  applied  by  Van  Dyke  demonstrates  the  importance  of  second  order  terms. 
The  metric  influence  of  curvature  has  been  given  by  Schultz-Grunow  arid  Brewer11  for 
laminar  and  incompressible  flows. 

Cebeci,  Hirsch,  and  Whitelaw12  analyzed  the  turbulent  boundary  layer  on  a  (convex) 
longitudinally  curved  surface.  In  their  treatment  they  replaced  the  laminar  viscosity  with  an 
effective  eddy  viscosity.  The  eddy  viscosity  was  specified  to  be  the  standard  Cebeci  & 
Smith13  formulation  modified  by  Bradshaw’s14  correction  for  longitudinal  curvature.  Eghlima 
et.al.15  used  a  similar  approach.  However  in  all  these  approaches,  the  assumption  has  been 
made  that  the  vorticity  vanishes  along  the  edge  of  the  boundary  layer,  i.e.,  the  outer 
flow  is  irrolational. 

It  has  been  found  that  the  standard  JANNAF  method2  for  predicting  boundary  layer 
losses  is  not  sufficiently  accurate  for  computing  boundary  layer  losses  for  space  engines. 
These  engines  exhibit  substantial  boundary  layer  growth  due  to  their  very  high  expansion 
ratio*.  For  small  engines,  a  major  portion  of  the  nozzle  flow  can  be  completely  enveloped  by 
the  wall  shear  layer  (up  to  90%*). 

The  two  major  deficiencies  of  the  JANNAF  boundary  layer  thrust  loss  calculation 
method2  are  that  the  effects  of  both  transverse  and  longitudinal  curvature  are  assumed  to  be 
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pmall.  Since  in  highly  expanded  flows,  the  boundary  layer  thickness  can  become  a  significant 
fraction  of  either  the  transverse  or  longitudinal  curvature,  this  effect  should  not  be  ignored. 

The  objective  of  this  paper  is  (1)  to  assess  the  significance  of  nozzle  wall  curvature 
pffects,  and  (2)  to  assess  performance  losses  due  to  viscous  effects  using  a  more  complete 
description  of  the  viscous  shear  layer  within  the  boundary  layer  equations  for  a  wide  range  of 
Space  applications  engines. 
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THEORY 


For  a  compressible  boundary  layer  flow  in  an  axisymmetric  nozzle,  the  governing 
equations  including  ‘he  transverse  and  longitudinal  curvature  effects  are  given  as12,15,17,18. 

Continuity: 


r 

(pur")  +  ^  pvF"  (1+ky) 


=  0 


(1) 


Stream  wise  Momentum: 


n  u _  8u  —  8u  pkuv  _  _  1 _  8p  1 

P  1+ky  9x  ^  3y  1+ky  1+ky  3x  +  r"(l+ky) 


8 

8y 


r"(l+ky) 


L^I-p73' 


;  pv  =  pv  +  p  v 


(2a) 


Normal  Momentum: 


ku2  _  1  8p 
1+ky  p  3y 


(2b) 


Energy: 


pu  9H  .  —  8H 
T+ky  83T  +  Pv  87  = 


1 


r"(l+ky)  *  3y 


r"(l+ky) 


[fe 


9H  +  u 
*  8y  +  V- 


1  - 


1  1 
P? 


du 


u%-pW7 


(3) 


The  influence  of  longitudinal  curvature  on  the  mean  boundary  layer  flow  is  modeled  by 
inclusion  of  the  centrifugal  force  term  in  the  nomial  momentum  equation9-11.  (See  i  quation 
2b). 


Figure  1  shows  the  definition  of  r,  r0,  R0  and  R. 
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Transverse  curvature  is  defined  as: 


_  y  cos  <}>  _  r_ 
r0  r0 

and  die  longitudinal  curvature  as: 
jj-  =  k(x) 


and 


h =  1  +  k(x)y 


(4a) 


(4b) 


(4c) 


The  boundary  conditions  for  equations  (1-3)  are: 


at 

O 

n 

U  -  0,  V  =  Vw  (x), 

(5a) 

and 

T  =  Tw  (x)  or  qw  (x)  is  known; 

(5b) 

at 

y  =  ymax» 

u  =  uc,  H  =  Hc,  p=pe  +  f^y  ulymax. 

(5c) 

The  boundary  condition  for  p(ymax)  was  taken  from  Van  Dyke9-10  who  matched  the 
"Inner  Limit  of  the  Outer  Expansion"  (j)lh  iteration  to  the  (j— l)*  iteration  of  the  "inner  flow" 
(i.e.,  the  boundary  layer).  Van  Dyke  showed  that  equations  (2a)  and  (5c)  include  all  second 
order  terms  caused  by  longitudinal  curvature  effect. 


To  include  the  longitudinal  curvature  effect  in  the  eddy  viscosity  model,  the 
Bradshaw’s19,20  expression  has  been  employed  to  correct  the  inner  eddy  viscosity  term  by 
multiplying  this  expression  by  S2  where: 


S  = 


T+K-]5 


(6) 


where  Rj  is  analogous  to  Richardson  number  and  p  is  reported  to  be  7  for  convex  and  4  for 
concave  surfaces19,20. 
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The  streamwise  momentum  equation,  (2a),  requires  8p/8x  along  lines  of  constant  y. 


Surfaces  of  constant  Tj  are  all  parallel  to  the  curved  wall. 
dt' 


"edge",  the  term 


I 


is  given  by: 


d£l 

X 


dur 


m 

m 

|3yJ 

X 

l^nJ 

X 

[8xJ 

At  the  boundary  layer 


(7) 


The  conventional  JANNAF  method2  for  evaluation  of  the  viscous  thrust  loss  in  rocket 
engine  nozzle  employs  both  the  boundary  layer  momentum  and  displacement  thicknesses. 
This  method  can  be  arrived  at  by  two  different  approaches  as  shown  by  Alber21.  The  basis  for 
the  approach  is  to  compare  the  thrust  of  an  inviscid  nozzle  to  a  viscous  nozzle  with  the  same 
mass  flow  rate.  In  the  equation  given  below,  the  first  term  represents  the  pressure  forces 
acting  on  an  inviscid  nozzle  and  the  second  and  third  terms  represent  the  total  stress  forces 
acting  on  the  viscous  nozzle. 

s  s 

AT  =  J  (p*— Peo)27trp  ^  dx  -  j*  (pw-pM)27t0y,5*cos<i>) .  (8) 

‘  0  0 


— Tpt^9?.S(}>)  dx  +  J  tw  27t(rp+5*cos<(>)cos<{)  dx 
o 


where  rp  is  the  distance  from  the  axis  to  the  potential  wall  and  p*  is  the  pressure  at  rp,  as  is 

C 

shown  in  Figure  1.  For  thin  boundary  layers,  the  conventional  assumption  is  to  let  rp=r0. 
However,  to  maintain  generality,  Equation  (8)  distinguishes  between  the  potential  and  real 
wall.  The  integral  form  of  the  momentum  equation  (wit*'  transverse  and  longitudinal  curvature 
effects)  can  be  derived  as: 


8 

8* 


+  PcUcTo 


pu  )  r_ 
pcUc.  r0 


dy 


~ro 


= 


(9a) 
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where:  +  PcUc^  =  e(x) 


(9b) 


The  details  of  the  derivation  can  be  found  in  Reference  16.  Alternate  forms  of  the  momentum 
and  displacement  thicknesses  for  axisymmetric  flows  are  defined  as: 


(momentum  thickness) 


_£U_ 

pcuc 


i_y_ 


dy 


(10a) 


(displacement  thickness) 


(10b) 


Implementing  Eqs.  (10)  in  Eq.  (9)  yields 


(ID 


Inserting  Eq.  (1 1)  in  Equation  (8)  and  neglecting  terms  of  order  52*,  the  thrust  deficit  is: 


J  (prpw]rp3x*dx“J  [pw"p«]  (rp8'0cos<t>)dx 


dr 


s 

+  /r0 
o 


-k(x)  f  puvdy  +  e(x)(5*-8+S2/2r0) 
J  r0 


.  cos<{>  dx 


I  COS(}). 


(pcuir05  -  r05*  ^  (pc  -  pj 


dx 


(12) 


Here  p*  is  the  pressure  at  the  edge  of  the  potential  flow  and  pc  is  the  pressure  at  the  edge  of 

C 

the  boundary  layer,  i.e.,  at  rw  -  8  cos  <J). 
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To  this  point,  the  thrust  loss  equation  has  been  derived  with  very  few  of  the  traditional 
boundary  layer  assumptions  being  madr  To  examine  the  effects  of  these  assumption^,  they' 
will  be  applied  one  at  a  time. 

1)  Non-dissipative  edge  condition,  i.e.,  e(x)*0.  In  the  absence  of  dissipative  or  applied 
forces  on  the  core  flow,  Bernoulli's  equation,  is  satisfied'  and  e(x)=0  (see'  e'q:  9b). 
Examples  of  flows  where  e(x)?K)  would  be  two  phase  core'  flow  of  MHD'-  flow.  The* 
second  term  of  the  third  integral  in  Eq.  (12)  disappears  when  e(x)=0'. 

2)  The  boundary  layer  is  thin  compared  to  the  local  radius  of  curvature;  he.,-  5/F?0«<l-. 
This  assumption  leads  to  the  following: 


3p/3y  =  0  or  p  =p  =p*=p 

wee 


and  g|  -  0 

;=  /  cos<t>-  [poUcr0^r05*(Pe-P^)']l  d* 

0 

-  /  (pc-pM)cos<})  (rp5*)  dx<  ('1 3) 


Integrating  the  first1  term1  in-  Eq.  (1*3)  by  parts  and-  assuming,  that  5?(s=o)  -  $(s=o)=0 
yields 


=  p  cu2rp§cos<H-p5’(p  0-p  Jcos<{> 


o  -y, 

0> 


In 

roJ 


'dx 


(1*)» 


3)  For  adiabatic  and  cooled  wall  flows,  the  displacement  thickness,  5*,  is  usually  less  than 
the  boundary  layer  thickness,  8.  The  assumption  applied  here'. is  that  5*c6s  <j)yr^<<r  1 
which  is  less  restrictive.  thaniS/r^«  1. 
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The  remaining  integral  in  eq.  (14)  vanishes  by  noting  that  _ 


5* 

rp/r0  =l-r-  cos^  =  1 

1  O 

To  the  same  order  of  accuracy,  Eq.  (14)  can  be  rewritten  as 

A  HP 

=  Pcucr0§cos<}>  -  r05*(pc-pJcos<{)  ( 1 5) 


4)  The  above  equation  reduces  to  the  standard  JANNAF2  relationship  for  thrust  loss  with 
the  additional  assumption  5cos<{>/r0«l.  That  is 

L  =  1  -  Scos<l)  ~  1 
r0  r0  “ 

which  leads  to 


9  =  &  and  8*  =  5* 


and  (15)  becomes: 


AT  =  2k  cos4> 


p  cu^0e  fo8 


(16) 


In  conclusion,  it  can  be  seen  that  in  order  to  reduce  the  general  form  of  the  thrust  loss 
Eq.  (8)  to  the  standard  JANNAF  form^,  i.e.,  Eq.  16,  requires  four  assumptions.  The  first 
assumption,  i.e.,  that  Bernoulli's  equation  is  satisfied  at  the  edge  of  the  boundary  layer  is  met 
in  most  liquid  rocket  engines  of  interest. 

The  other  three  assumptions  deal  with  varying  degrees  of  how  thin  the  boundary  layer 
is  compared  to  a  given  length  scale.  The  "thin"  assumptions  are: 

0  The  boundary  layer  is  thin  compared  to  the  local  longitudinal  radius  of 
curvature,  i.e.,  8/R0  «  1. 
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o  The  displacement  between  the  viscous  and  potential  wails  is  small  compared  to 
the  local  wall  radius,  i.e.,  5*cos<j>/r0  «  1,  which  leads  to  rp  ~  r0. 

o  The  boundary  layer  thickness  is  small  compared  to  the  local  wall  radius,  i.e., 
6cos<j)/r0  «  1,  which  leads  to  r  ~  rp  ~  r0. 

RESULTS  AND  DISCUSSION 

The  BLM  (Boundary  Layer  Module)  within  the  TDK  (Two-Dimensional1  Kinetics) 
code  was  modified  to  include  the  normal  momentum  equation  together  witir  longitudinal' 
curvature  terms. 

Solutions  to  the  momentum  and  energy  equations-  can-  be  obtained*  in*  a  veiy  efficient 
manner  by  using  the  block  elimination  method  as  discussed-  By  Kelleri  The  iriiplicit  finite 
difference  scheme  that  has  been-  developed  by  Keller  and*  Cebeci22  was  implemented'  to  obtain' 
numerical  solution  to*  these  equations16'.  The  BLM;  Code  was  Validated  By  comparison  ta 
another  boundary  layer  code,.  MABL23.  The  code  was  also  verified*  for  flow  over  a  flat  plate 
and  available  experimental  data^  for  flow  over  a  convex  surface.  These  results  are  presented  in 
Ref.  16. 

The  nozzle  wall  is  not- generally  known  analytically,  and  is  often  obtained  from  tabular 
input  data  by  spline  fitting.  Calculation  of  longitudinal  curvature  which  requires  the  second 
derivative  of  a  nozzle  wall  creates  "severe"  fluctuations  of  the  second*deriVative  which  in  turn 
can  cause  numerical  instabilities;  It  was  decided  to  divide  the  nozzle  wall  info  numbers  Of 
equally  spaced  intervals  and  let  a  circle  pass  through  every  three  points.  The  radius  of  the 
circle  is  taken  to  be  the  radius  of  the  curvature  at  the  middle  point.  This  method  has  been 
verified  against  an  analytical  wall  and  the  results  are  satisfactory16 ,24. 

The  importance  of  longitudinal  curvature  is  determined  by.  the  magnitude  of  5/R0,  i.e:, 
the  latio  of  the  boundary  layer  thickness  to  the  radius  of  curvature.  Figs.  2a  arid  f2B' indicate 
this  value  for  ASE  and  SSME  nozzles,  respectively.  The  magnitude  of  8/R0  for  SSME  is* 
larger  than  ASE  nozzle  and  thus  a  larger  pressure  gradient  across  the  boundary  layer  is 
expected.  This  can  be  observed:  from  a  comparison  of  velocity  profiles  in  Figs.  3  arid  4. 

The  effect  of  longitudinal  curvature  on  performance  is  depicted  in  Table  1.  It  can  be 
seen  that  although  the  longitudinal  curvature  affects  the  mean-  flow  velocity,  temperature  and ' 
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pressure  profiles  in  the  boundary  layer,  it  has  insignificant  impact  on  nozzle  performance. 
However,  the  methodology  being  used  to  evaluate  viscous  loss  has  to  distinguish  oetween  a 
thin  and  thick  boundary  layer. 

Essentially,  Eq.  (A)  in  this  table  which  is  the  standard  JANNAF  method  is  not  valid  for 
thick  boundary  layers,  because  of  the  magnitude  of  term  r/r0  where  r  is  the  distance  from  the 
axis  of  symmetry  to  the  edge  of  the  boundary  layer.  As  the  value  of  r/r0  increases,  the  thin 
shear  layer  assumption  yields  erroneous  results  for  boundary  layer  losses,  i.e.,  Eq.  (A).  Eq.  (B) 
was  developed  for  axisymmetric  flows,  and  for  the  thick  boundary  layers  is  expected  to  yield 
quite  different  values  from  Eq.  (A)  and  it  can  be  seen  that  for  nozzles  with  high  area  ratio  such 
as  ASE,  RL-10  and  XLR134,  the  difference  between  these  equations  is  more  significant.  This 
difference  is  due  to  a  thick  boundary  layer  when  r/r0  is  significantly  less  than  unity.  To  be 
more  exact  rp  should  be  replaced  for  r0  in  Eq.  (B)  as  shown  in  the  analysis  and  upon  doing 
that  in  fact  the  effects  are  more  severe  on  the  high  area  ratio  nozzles,  Eq.  (C).  The  lesults  of 
adding  another  higher  order  term  to  Eq.(C)  has  been  indicated  in  Eq.  (D).  It  can  be  seen  that 
addition  of  this  term  makes  a  change  of  up  to  about  8%  in  the  thrust  loss  for  high  area  ratio 
nozzles.  Using  the  integrated  wall  shear  method  i.e.,  Eq.  (E)  yields  results  quite  different  (up 
to  16%)  from  Eq.  (A)  for  ASE  nozzle,  (area  ratio=  400),  RL-10  (area  ratio  205)  and  XLR134, 
(area  ratio  767:1).  However,  as  the  higher  order  terms  are  added  to  Equation  (D),  the  results 
become  closer  to  values  of  Eq.  (E).  This  is  expected  because  of  mathematical  equivalency  of 
both  equations.  The  overall  results  indicate  that  the  standard  JANNAF  method  for 
performance  prediction  is  not  adequate  and  can  produce  fairly  large  error  performance 
predictions  for  high  expansion  nozzles. 

The  results  indicated  in  this  table  justifies  the  usage  of  Eq.  (E)  for  the  thrust 
Calculations  and  Eq.  (E)  should  be  used  to  obtain  the  thrust  loss.  Furthermore,  the  inclusion 
of  longitudinal  curvature  is  not  crucial  to  performance  prediction  in  contoured  wall  of  a  rocket 
engine  nozzle. 


CONCLUDING  REMARKS 

The  results  in  this  work  indicated  that  although  the  effects  of  the  longitudinal  curvature 
are  of  importance  on  the  mean  flow  velocity,  pressure  and  temperature  profiles  in  the 
boundary  layer  in  strongly  curved  nozzle  contours,  it  does  not  have  a  significant  impact  on  the 
overall  prediction  of  nozzle  performance.  However,  for  nozzles  with  a  thick  viscous  layer,  the 
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standard  JANNAF  method  which  is  formulated  for  thin  shear  layers  "over-estimates"  the 
boundary  layer  thrust  deficit  and  it  should  be  replaced  by  the  equation  (E)  in  Table  1  which 
will  be  referred  to  as  "wall  shear  method". 
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Table  1  Boundary  Layer  Thrust  Deficit  for  Variety  of  Engines  in  the  Absence 
and  Presence  of  the  Longitudinal  Curvature.  Cold  Wall  Case 
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NOZZLE 

THRUST  DEFICIT,  lbf.  aoc 

./  ltxn 

WITHOUT  LONGITUDINAL  CURVATURE 

WITH  LONGITUDINAL  CURVATURE 

EQ.  A  EQ.  D 

EQ.  C 

EQ.  D 

EQ.  E 

EQ.  A 

EQ.  8 

EQ.  C 

EQ.  D 

EQ.  E 

SS.  E 
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6.3184 

6.4758 

6.0387 
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ASE 

14.1985  13.6977 
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BOUNDARY  LAYER  LOSS  MODELS  IN  NOZZLE  IN  PERFORMANCE  PREDICTION 
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Software  and  Engineering  Associates,  Inc. 
Carson  City,  Nevada 


ABSTRACT 

This  paper  details  the  importance  of  boundary  layer  shear  stress  models  in  predicting  the  boundary  layer 
parameters  of  interest.  Laminar  and  turbulent  boundary  layers  arc  investigated.  Different  shear  stress  models  arc 
considered  for  a  wide  range  of  liquid  and  solid  propellant  systems  and  over  a  range  of  nozzle  expansion  ratios.  Simple 
empirical  fits  to  the  computed  results  models  for  solid  and  liquid  systems  arc  also  given. 


INTRODUCTION 


The  trend  toward  very  large  area  ratio  nozzles,  which  result  in  performance  gains  for  space  propulsion 
applications,  has  increased  the  need  for  detailed  knowledge  of  the  momentum  losses  due  to  nozzle  viscous  effects  (i  c  , 
boundary  layer).  These  losses  degrade  overall  system  performance,  such  as  increasing  system  weight,  decreasing  useful 
payload  weight,  and/or  decreasing  effective  system  range.  Another  important  factor  in  the  designing  of  propulsive 
nozzles  is  the  detailed  knowledge  of  heat  transfer  at  the  wall  for  rcgencrativcly  cooled  walls  and/or  material  performance. 

Because  of  the  importance  to  rocket  propulsion,  the  Astronautics  Laboratory  at  Edwards  Air  Force  Base  has 
sponsored  the  Boundary  Layer  Study  Contract  to  improve  the  understanding  and  computational  predictive  capabilities  for 
boundary  layers  in  rocket  nozzles  with  high  area  ratios.  As  part  of  this  study,  the  sensitivity  of  the  computed  boundary 
layer  loss  to  various  parameters  was  investigated.  The  effects  of  both  longitudinal  and  transverse  curvature  were  reported 
in  Reference  1  and  2.  Interactions  between  the  viscous  wall  shear  layer  and  the  inviscid  core  have  been  reported  in 
Reference  3.  This  paper  discusses  the  importance  of  validating  the  turbulent  shear  stress  model. 


The  methods  used  to  calculate  the  boundary  layer  performance  loss  are  presented  first.  The  magnitude  of  the 

differences  between  methods  of  computing  Al  are  not  small2,4  and  should  be  addressed  by  the  propulsion 

spBL 


community. 


Results  of  boundary  layer  losses  arc  then  presented  for  a  variety  of  liquid  and  solid  systems.  Computations  arc  for 
both  laminar  and  turbulent  flows  (eddy  viscosity  model)  using  the  TDK5,  SPP6,  and  VIPER7^  codes.  Results  comparing 

eddy  viscosity  and  k-€  turbulence  models  arc  presented  for  one  case. 


NOMENCLATURE 


I  Specific  impulse 

K  Kinetic  energy  of  turbulence 

L  Unit  length 

p  Pressure 

Pr  Prandti  number 

R^,  Radius  of  longitudinal  curvature 

r  Radius  of  transverse  curvature 

Re  Reynolds  number 

t  =  y  cos<{>/r0 

T  Thrust 

u,v  Velocity  components 


CreckXettcrs 

6  Boundary  layer  thickness 

5#  Displacement  thickness 

e  Turbulence  dissipation  rate 

also  expansion  ratio 
<{>  Wall  angle 

p  Viscos.cy 

v  Kinematic  viscosity 

p  Density 

8  Momentum  thickness 

x  Shear  stress 


SubsciiEis 

BL  Boundary  layer 

D  Delivered 
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CALCULATION  OF  THE  BOUNDARY  LAYER  LOSS  IN  NOZZLE  FLOWS 

The  conventional  JANNAF  method  9~10  for  evaluation  of  the  viscous  thrust  loss  in  rocket  engine  nozzle  employs 
both  the  boundary  layer  momentum  and  displacement  thickness.  This  method  can  be  arrived  at  by  two  different 
approaches  as  shown  by  Albcr11.  The  basis  for  the  approach  is  to  compare  the  thrust  of  an  inviscid  nozzle  to  a  viscous 
nozzle  with  the  same  mass  flow  rate.  In  the  equation  given  below,  the  first  term  represents  the  pressure  forces  acting  on 
an  inviscid  nozzle  and  the  second  and  third  terms  represent  the  total  stress  forces  acting  on  the  viscous  nozzle. 

AT  n  1  (p*  -  pJ2*rp  dx  - 1  (pw  -  PJ2n(rp+5*cos<j»  dx  (1) 

0  c  0 


+  J  tw  2nfrp+6*cos<j))cos<j)  dx 
o 


where  rp  is  the  distance  from  the  axis  to  the  potential  wall  and  p*  is  the  pressure  at  rpt  as  is  shown  in  Figure  1. 

When  the  boundary  layer  characteristics  lengths,  8,  0,  and  5*,  are  assumed  to  be  small  compared  to  the 
characteristic  nozzle  lengths,  r  and  R0  then  equation  (1)  reduces  to  the  standard  JANNAF  relation10 


AT  =  2k  cost})  pcuJro9  -  r06*(pc  -  p^) 


(2) 


If  the  boundary  layer  thickness  assumption  that  8/Rc  «  1  is  relaxed,  then  the  modified  JANNAF  relation  can  be 
derived4,  i.c., 


=*  peu|ru0cos4)  -  roS*(pe  -  p°*)co  sty  (3) 

using  the  alternate  forms  of  the  momentum  and  displacement  thicknesses  for  axisymmctric  flows  shown  below 

(momentur.  thickness)  0  =  f  ~  fl  -  dy 

ro  Pcuc  l  uej 


(4a) 


j 

(displacement  thickness)  5*  =  J  —  1  -  dy 

or0  PeUeJ 


(4b) 


Table  1  (from  Reference  2)  shows  that  the  differences  in  boundary  layer  loss  as  computed  from  Equation  1,  2,  and 
3  (columns  E,  A,  and  B  respectively)  are  significant  for  nozzles  with  thick  boundary  layers. 


Figure  1.  System  Geometry 


X 


Table  1  Boundary  Layer  Thrust  Deficit  for  Variety  of  Engines  In  the  Absence 
and  Frasanca  of  the  Longitudinal  Curvatura.  Cold  Wall  Cate 


THRUST  DEFICIT,  lbf.  sec./ltxn 


KOZZLE  WITHOUT  LONGITUDINAL  CURVATURE  WITH  LOHGITUDIHAL  CURVATURE 


EQ.  A 

EQ.  B 

EQ.  C 

EQ.  D 

EQ.  E 

EQ.  A 

EQ.  B 

EQ.  C 

EQ.  D 

EQ.  E 

SSME 

6.6843 

6.7730 

6.7650 

6.2088 

6.2233 

6.3720 

6.3184 

6.4756 

6.0387 

6.5036 

ASE 

14.1985 

13.6977 

13.4572 

12.3531 

11.9394 

14.2351 

13.7440 

13.5000 

12.3782 

11.7659 

RL-10 

11.6549 

11.0671 

10.6094 

10.4630 

10.4930 

11.7732 

11.0714 

10.8591 

10.4517 

10.4964 

XDELTA 

1.0459 

1.8160 

1.7902 

1.7220 

1.6132 

1.6450 

1.8160 

1.7670 

1.7215 

1.8211 

BC4515 

0.0688 

0.8534 

0.8481 

0.8142 

0.8373 

0.6686 

0.8534 

0.8221 

0.8141 

0.8378 

BC1010 

1.8578 

1.7941 

1.7866 

1.8000 

1.8790 

1.8582 

1.7942 

1.6121 

1.8001 

1.8788 

IUS 

1.0064 

1.9523 

1.9223 

1.8381 

2.0943 

1.9862 

1.9521 

1.8524 

1.8360 

2.0628 

XLR-134  27,3262 

25.6350 

25.3452 

24.6576 

23.2534 

27,3525 

25.6539 

25.2531 

24.6322 

23.1904 
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in  all  tha  above  equations  T  (x)  «  0 
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ROLE  OF  TURBULENCE  MODELS 


The  unsteady  Navicr-Stokcs  equations  arc  generally  recognized  as  being  capable  of  accurately  describing  both 
laminar  and  turbulent  flows.  However,  the  difficulties  associated  with  resolving  the  various  length  scales  in  turbulent 
flows  has  removed  direct  solution  of  these  equations  from  practical  engineering  usage.  Instead,  the  time  or  Reynold’s 
averaged  Navier-Stokes  equations  are  solved  for  cases  of  engineering  interest.  The  time  averaging  of  the  equation 
generates  "apparent”  stresses  and  heat  fluxes  which  are  modeled  empirically.  This  empiricism  relates  mean  flow 
variables  to  the  apparent  stresses. 

Two  classes  of  turbulence  modeling  are  currently  popular.  The  first  class  is  the  eddy  viscosity  or  mixing  length 
methods  which  trace  their  origins  to  Boussinesq  and  Prandtl.  These  methods  arc  also  referred  to  as  simple  algebraic  or 
zero  equation  models.  The  second  class  uses  transport  equations  to  evaluate  the  turbulent  Reynolds  stress  terms.  A 

frequently  used  twoequation  model  is  the  k-e  method  which  was  first  proposed  by  Harlow  and  Nakayama12.  Some 

popular  k-€  methods  follow  the  work  of  Jones  and  Launder12  and  Launder  and  Spalding14. 

The  JANNAF  standard  codes  for  evaluating  boundary  layer  losses,  i.c.,  SPP-BLM6,  TDK-BLM/MABL5,  and 
BLIMP-J15  all  use  the  Cebeci-Smith  eddy  viscosity  model16.  The  VIPER  code  being  developed  under  this  effort  includes 

both  the  Cebeci-Smith  eddy  viscosity  model  and  a  k-e  model,  see  Reference  7  and  8. 

With  the  success  of  viscous  flow  calculations  it  is  easy  to  forget  the  origin  and  amount  of  empiricisms  that  go  into 
these  calculations.  All  of  these  turbulence  models  were  developed  for  air  and  pressure  gradients  which  seem  very  modest 
by  rocket  nozzle  standards.  A  study  conducted  by  Evans1'  comparing  3  similar  eddy  viscosity  models  (Kendall, 
Bushnell,  and  Cebeci-Smith)  concluded  that  large  differences  in  heat  transfer  could  occur  between  these  models.  The 
main  culprit  was  assumed  to  be  the  way  which  these  models  handle  density  and  temperature  gradients  across  the 
boundary  layer,  especially  in  the  law  of  the  wall  region. 

Basclinc.CQmputaiiQns 

To  establish  a  baseline  for  the  computations,  a  series  of  four  liquid  engines  (RL10,  SSME,  XLR134,  and  ASE) 
were  run  on  the  TDK/BLM  codes  and  eight  solid  motor  eases  on  SPP/BLM.  The  results  of  these  computations  are  shown 
in  Figure  2.  The  major  reason  for  higher  losses  on  liquid  systems  is  that  these  systems  tend  to  have  cooled  walls  and 
hence  higher  boundary  losses.  No  attempt  was  made  to  correct  for  the  effects  of  regen  heat  addition  to  the  core  flow  on 
the  overall  boundary  layer  loss.  As  a  result  of  these  computations,  we  found  the  following  expressions  to  be  reasonable 
fits  of  the  data. 


solid  propellant  systems 

AT  /I  x  100  «  0.32  +  3.887  x  HHe 
S?BL  SPD 

liquid  propellant  systems 

AT  /I  x  100  =  1.065  +  6.192  x  10*  e 

spBL  SPD 

where  t  is  the  nozzle  area  ratio 

While  the  above  can  be  used  for  estimates  of  the  boundary  layer  loss,  care  should  be  exercised  for  liquid  systems 
since  the  amount  of  heat  extracted  from  the  flow  (and  hence  AI  )  varies  greatly  from  engine  to  engine  and  is  not  a 

spBL 

function  of  area  ratio. 

ENSJtaults 


Computations  were  performed  for  six  !iqi»;d  propellent  rocket  engines  using  ihe  VIPER7-#  code.  The  VIPER  code 
was  selected  for  this  study  because  it  contains  all  of  the  models  to  do  a  complete  and  insistent  set  of  calculations.  It  can 

treat  both  laminar  and  turbulent  flow  with  reacting  chemistry.  It  has  tw  i  turbulence  models,  Cebeci-Smith  and  k-e,  and 
can  handle  very  thick  boundary  layers  since  the  core  flow  and  wall  shear  layer  are  directly  coupled.  The  boundary  layer 
loss  was  computed  as  the  integral  of  the  shear  stress  along  the  wall  (last  term  in  equation  1). 

Laminar  and  turbulent  (eddy  viscosity)  results  were  obtained  for  all  six  nozzles.  The  characteristics  of  these 
engines  are  given  in  Table  2  and  the  results  are  shown  in  Figure  3.  While  there  is  a  significant  amount  of  scatter  in  the 
computed  results,  it  is  clear  that  the  turbulent  losses  arc  from  2  to  4  times  larger  than  the  laminar  losses. 


6.00 


O  U0UI0  STSTErtS 


D  SOHO  5T5TEOS 


Figure  2.  Boundary  Layer  Lots  as  a  Function  of  c  Area  Ratio 


Tabic  2.  Engine  Characteristics 


Engine  Name 

Chamber 

Pressure 

Throat 

Radius(in) 

Expansion 

Ratio 

NASA/Lewis  Hi-E 

360 

.5 

1025. 

XLR-134 

510 

.396 

767.9 

RCS 

150 

1.021 

28.46 

SSME 

3285 

5.1527 

77.5 

RL  10 

394.3 

2.57 

205.03 

ASE 

228 

1.254 

400.7 

Unfortunately  the  k-€  turbulence  model  incorporated  in  VIPER  turned  out  not  to  be  as  robust  as  we  had  hoped  and 
we  were  able  to  obtain  reasonable  results  for  only  one  case.  That  case  was  the  XLR-134  high  expansion  ratio  nozzle. 
The  following  table  shows  the  results  for  this  nozzle  using  frozen  chemistry. 


Tabic  3.  XLR  134  Boundary  Layer  LOSS  Results 


Turbulence  Model 

Delivered 

!sp 

A1 

sPbl 

VsP% 

Cebeci  Smith 

456.377 

17.839 

3.909 

k-€ 

452.833 

23.688 

5.231 
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No  matter  how  the  comparison  of  the  above  is  made,  a  one  percent  difference  in  delivered  I  has  to  be  considered 

sp 

significant. 

The  magnitude  of  the  boundary  layer  loss  and  the  discrepancy  between  turbulence  models  can  only  be  resolved  by 
very  careful  boundary  layer  measurements  in  an  environment  similar  to  a  rocket  engine. 

It  is  obvious  from  the  results  presented  here  that  experimental  validation  of  boundary  layer  loss  models  is  required 

if  the  design  margins  on  future  high  area  ratio  space  systems  arc  to  be  lowered.  To  put  this  statement  in  perspective,  the 

rule  of  thumb  is  that  a  loss  of  1%  in  I  results  in  a  2%  loss  in  payload. 

sp 

CONCLUSIONS 


The  magnitude  of  the  boundary  layer  performance  loss  has  been  shewn  to  range  1  to  6%  of  the  total  delivered  ISp 

for  liquid  propellant  rocket  engines  with  moderate  to  high  area  ratios.  The  effect  of  turbulence  has  been  shown  to 
increase  the  performance  loss  over  a  laminar  boundary  layer  by  factors  of  2  to  6  making  the  onset  of  turbulence  a  critical 
factor  for  rniall  high  area  ratio  nozzles. 

The  uncertainty  in  turbulence  models  had  been  previously  estimated  to  be  in  the  10  -  25%  range  on  predicted  I 
loss.  The  calculations  presented  here  showed  u  range  uf  25  -  33%  depending  on  the  *c*vd  base. 
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ABSTRACT  I.  INTRODUCTION 


The  effects  of  thick  boundary  layers  on 
nozzle  flow  performance  has  been  Investigated. 

The  relationship  between  the  boundary  layer  and 
the  core  flow  for  a  variety  of  real  engine 
parameters  has  been  established.  Eulcr/boundary 
layer  solutions  are  compared  against  the  full 
Navler-Stokes  solver  and  PNS  solver.  These 
solutions  have  been  compared  to  the  standard 
JANNAF  TDK/BLM  computer  code. 

The  interaction  between  the  viscous  layer  and 
the  core  flow  has  beon  studied  for  a  thick 
boundary  layer  where  the  '’thin*'  shear  layer 
equations  are  not  valid.  This  study  Is  for 
perfect  gas  with  constant  y. 

The  validity  and  restrictions  for  "thlnN 
shear  layer  assumption  for  thick  boundary  layers 
has  been  studied. 


ABBREVIATIONS 


ASE  Rocketdyne  Advanced  Space  Engine 

BUI  Boundary  Layer  Module  of  the  TDK/BLM 

Computer  Program 

CFL  ~ourant,  Levy,  and  Friedricks  (Number) 

I  Specific  Impulse 

OTV  Orbital  Transfer  Vehicle 

SBS-1A  Hughes  Spacecraft  Engine 

SSHE  Space  Shuttle  Haln  Engine 

TDK  Two-Dimensional  Kinetic  Computer  Code 

XLR-134  Low  Thrust  Cryogenic  Engine,  OTV 
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The  trend  towards  very  large  area  ratio 
nozzles,  which  result  In  performance  gains  for 
space  applications,  has  Increased  the  need  for 
detailed  knowledge  of  the  momentum  losses  due  to 
nozzle  viscous  effects  (i.e.,  boundary  layer)  and 
the  flovfleld  parameters  in  propulsion  systems. 
These  losses  degrade  overall  system  performance 
•s  it  directly  affects  payload  efficiencies  such 
as  increasing  system  weight,  decreasing  useful 
payload  weight,  and  decreasing  effective  system 
range. 

As  the  interest  In  unders*  tiding  of  high 
speed  flows  and  performance  of  propulsive  nozzles 
Increases,  there  Is  need  for  the  development  of 
more  precise  computational  techniques.  The 
performance  loss  .in  high  expansion  ratio  nozzles 
due  to  the  viscous  layer  is  m  critical  design 
parameter  In  t*ic  nozzle  geometry. 

Recent  studies'  '  Indicate  that  for  high  area 
ratio  nozzles,  the  viscous  layer  can  become  quite 
thick  such  that-  the  evaluation  of  losses  due  to 
this  layer  basea  upon  the  traditional  Prandtl 
thin  sheer  layer  approximation  could  be 
Insufficient  for  its  application.  Although  the 
second  order  boundary  layer  equations  could  be 
utilized  for  obtaining  improved  results,  there  Is 
still  a  need  for  superior  techniques  to  treat 
"thick"  viscous  layers. 

Very  powerful  computational  capabilities  have 
been  developed  in  the  application  of  finite 
difference  techniques  to  the  solution  of  the 
fluid  dynamic  equations  governing  compressible 
flows . 

The  most  notable  algorithm  was  introduced  by 
MacCormack  in  his  classic  1969  paper  on  the 
effect  of.ylscoslty  on  hypervelocity  Impact 
cratering''5  .  His  explicit  solution  technique, 
which  Is  a  variation  on  the  Lax-Ucndroff  second 
order  mothod,  uses  a  predictor -corrector  operator 
to  achieve  second  order  accuracy  while  using  a 
first  order  accurate  finite  difference  approach. 
Since  this  method  Is  easily  applied  to  complex 
flowfleld  problems ,  It  Is  widely  used  in  the  CFD 
community  today  . 


To  resolve  the  viscous  layer  portion  of  the 
flow  in  Navler-Stokes  solvers,  a  fine  mesh  should 
be  considered  near  the  wall  region.  This  will 
result  In  increasing  the  computation  time 
considerably.  To  some  extant,  this  difficulty 
has  been  overcome  by  the  development  of  FNS 
(parabolized  Navler-Stokes)  equations  and 
PNS-based  computer  programs. 


Tlie  parabol  izcd  Navler-Stokos  equations 
neglect  llie  strenmwlse  diffusion  term,  which 
along  with  special  treatment  of  the  subsonic 
region  of  the  boundary  layer  removes  the  spatial 
elliptic  I ty  from  the  steady  form  of  the  equations 
and  permits  a  solution  by  streainwlse  marching 
computation'll  techniques.  The  major  advantage  of 
•he  PNS  co^es,  in  general,  is  that  they  solve  a 
steady  form  of  the  governing  equations  and  are 
much  more  efficient  than  time ‘dependent  codes  for 
solution  of  the  unsteady  Nnvler- Stokes  equations. 
Due  to  these  attributes,  I’NS  codes  are  widely 
used  today  to  analyze  fluid  dynamics  of 
compressible  flow' 

Our  goal  In  this  work  was  to  establish  a 
fundamental  understanding  of  thick  viscous 
layers.  Supersonic  flows  in  propulsive  nozzles 
have  been  considered  ir*  this  work.  Three  codes 
(NS,  TNS,  and  boundary  layer  solvers)  were 
selected  and  exercised  on  the  candidate  engines 
embracing  laminar  to  fully  turbulent  boundary 
layer  flows. 


II.  TEST  CASES  AND  CODES 


A  list  of  candidate  engines  and  their 
specif Icatlons  Is  shown  In  Table  l.  The  test 
cases  cover  the  whole  range  of  laminar  (300:1 
nozzle)  to  completely  turbulent  nozzle  flows  (ASE 
and  SSME) .  The  software  selected  to  perforn0this 
study  were  the  TDK/BLM  codeU\  a  PNS  codeU  \ 
the  A3fj2DS  (Euler)  codeU  \  and  the  VNAP2  (NS) 
code'  . 

Some  of  the  computer  programs  had  to  be 
modified  to  Increase  the  resolution  and/or 
computation  accuracy.  Both  the  Invlscld  and 
viscous  solutions  are  obtained  and  compared. 


III.  1NVISCTD  SOLUTION 


A  comparison  between  the  TDK,  VNAP2 ,  and 
AXI2DS  code  results  for  the  nozzle  mass  flow 
rate,  thrust,  and  I  are  presented  in  Table  2. 
Tho  VNAP2  Code  (NS  IBlvor)  vns  modified  Co 
separate  the  subsonic/transonic  and  supersonic 
portion  of  the  flow.  The  validity  of 
“segmentation”  of  the  nozzle  lias  been  discussed 
elsewhere'  '  and  it  has  been  shown  that  this 
method  results  In  accurate  solutions  In  a  fairly 
reasonable  amount  of  computation  time. 

There  is  excellent  agreement  between  VNAP2  in 
the  Euler  mode  and  TDK.  AXI2DS  yields  comparable 
results;  however,  due  to  the  poor  grid  resolution 
in  AXI2DS,  the  results  are  slightly  different 
from  the  other  two.  From  this  table  |.t  can  be 
concluded  that  for  low  thrust  engines,  such  as 
the  Hughes  300:1  engine,  I  is  very  sensitive  to 
the  nozzle  mass  flow  rate  §Be  to  the  small  mass 
flow  rates. 

So,  for  low  thrust  engines,  It  is  necessary 
that  the  methodology  being  used  must  conserve 
mass  as  much  as  possible.  Also,  from  Table  2,  it 
can  be  observed  that  poor  grids  have  significant 
effect  on  the  higher  area  ratio  nozzles,  l.e., 
for  high  expansion  ratio  nozzles,  such  as  the 
XLR-13A,  an  enormous  number  of  grid  points  is 
necessary  to  solve  the  flowfteld  corrcctlv. 


IV.  VISCOUS  SOLUTIONS 


The  results  of  the  thrust  and  I 
calculations  for  all  the  test  cases^Snd  for  the 
viscous  flows  are  tabulated  in  Table  3.  The 
Initial  line  properties  for  the  PNS  Code  were 
taken  from  the  VNAP2  solution.  For  the  low 
thrust  engine,  the  flow  Is  laminar  and  the 
viscous  layer  occupies  about  301  of  the  nozzle 
(S  find  disagreement  between 

TDK/BiM  and  VNAP2  Is  basically  due  to  the 


Tabic  1:  Candidate  Engines  for  Thick  Boundary  Layer  Study 


ENGINE 

APPLICATION 

PROPELLANT 

TiffllCT 

<lbf) 

OIAMJEH 

MissufE 

(pala) 

THHOVT 
RADIUS,  R* 
(Indira) 

lOJCLD’S 

NllDUl 

nCH« 

EWANS  ICN 
RATIO 

S5M:  (Spooe  Shuttle 
Main  Engine) 

SST0 

Lox/ai2 

*163.000 

2935.7 

5.1527 

1.18  x  107 

77.5:1 

ASE  (Rocketdjne, 
Advanced  Space 

Engine) 

OIV 

uox/ai2 

22,600 

2287 

1.25*1 

2.20  X  I06 

*00:1 

ow’ia 

Spacecraft  Engine) 

Spacecraft 

¥n 

5.n’5 

10C> 

.0925 

3.50  x  iO4 

300:1 

XLfi-134 

orv 

lOX/Oi2 

511 

510 

.396 

1.80  x  I05 

761:1 
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TabU  2?  Comparison  of  Thrust  Calculations  for  TOC,  WAP2  and  AXI20S  for  on  tnvlscid  flou  of  a  Rarfect  Cat,  Adiabatic  Wall 


TDK 

VXAP2  (Cuter  Mode) 

AX 1 20 S  (Euler  .Mode) 

HASS 

MIDI 

MASS 

X  ISP 

MIDS 

70  (AXIAL)  X  44  (RADIAL) 

3USP 

tome 

new 

•ATE, 

THRUST , 
Ibf 

UP. 

icc.tbf 

AXIAL  x 
IA0IAL 

CF1 

flow 

RATE 

THRUST, 

Ibf 

ISP. 
tec. Ibf 

DIFFERENCE 

WITH 

HASS 

FLOW 

THRUST, 

Ibf 

IspT 

tec. Ibf 

DIFFERENCE 

WITH 

m 

Ibm/SfC 

-res— 

fURSOilTr" 

SUPERSONIC 

TRANSONIC 

SUPERSONIC 

I  be/ tec 

lb* 

TDK 

ttME, 

Un/sec 

lb« 

TDK 

■ 

9 

450764.6 

403.36 

5.22 

ESHC 

1051.95 

450690.9 

■ 

1059.96 

425.26 

.07 

1059.96 

H 

9 

22484.48 

456.90 

22195.73 

451.03 

1.23 

*SE 

49.35 

22535. 1) 

u 

151  x  51 

9 

49.2101 

-.05 

i9.2J 

iEH 

9 

5.441 

246.00 

.02212 

5.346 

241.67 

2.85 

gsjg 

0.02193 

5.455 

248.76 

n 

.02212 

1.11 

| 

1.13193 

520.94 

460.12 

1.13193 

481.03 

424.96 

7,40 

XLR-134 

1.11455 

511.49 

1 

-.28 

resolution  of  the  viscous  layer,  particularly 
near  the  wall.  These  results  ore  obtained  for 
variable  grid  spacing  with  more  gridlines  close 
to  the  wall.  However,  in  TDK/BLM  about  120 
points  are  assumed  in  the  viscous  layer  as 
opposed  to  about  20  points  in  the  NS  or  PNS 
solvers.  So,  It  can  be  concluded  that  the 
discrepancy  between  these  results  is  due  to  the 
resolution  of  the  viscous  layer. 


There  are  still  fundamental  questions  th»C 
have  to  be  answered:  that  is  how  thick  is  the 
subsonic  layer  and  how  accurate  can  an  NS  or  PNS 
solver  predict  the  thickness  of  this  sublayer? 

Two  of  the  candidate  engines,  ASE  and  XLR-134, 
were  chosen  to  answer  these  questions  for  further 
studies. 

In  ASE,  the  viscous  layer  occupies  about  15% 
of  the  nozzle  and  in  the  XLR-134  nozzle,  this 
value  is  about  20%  (£  qgjA)*  Both  adiabatic 
and  cooled  walls  are  considered;  and  subsonic 
layer  thickness,  together  with  a  number  of  points 
across  the  viscous  layer  and  subsonic  layer  are 
depicted  in  Table  4.  From  this  table,  it  con  be 
observed  that  generally,  If  an  NS  solver  is 
applied  to  supersonic  nr-^le  flow  to  resolve  the 
viscous  layer,  an  enormous  number  of  points  will 
be  needed  which  makes  the  computational  time 
prohibitively  large.  This  Cable  also  indicates 
that  the  subsonic  layer  thickness  is  about  1%  of 
the  boundary  layer  and  the  mass  travelling  in 
this  sublayer  is  so  small  such  that  It  cannot 
have  any  major  effect  on  the  performance.  The 
'thickness  of  the  subsonic  layer  Is  very  small 
compared  to  the  viscous  layer  and  the  flowfleld, 
and  does  not  establish  a  concern  about  the 
Information  travelling  backward  in  the  nozzle  In 
this  sublayer  due  to  the  elliptic  nature  of  the 
governing  equations. 


v.  shear  layer  core  flow  interaction 


In  the  nozzle  expansion,  at  some  point,  the 
Reynold's  number  can  become  low  enough  that  the 
thin  shear  layer  (*!$L)  or  boundary  layer 
equations  no  longer  yield  acceptable  accuracy  for 
performance  prvdicricn.  There  is  no  simple  way 
to  detect  when  the  75L  equations  are  no  longer 
adequate. 

There  are  several  studies  in  the  area  of  the 
invlscid-vlscous  interaction.  Among  these^are. 
the  methods  of  Carter^  '  \  U  Bflleur 

Wigton^  ,  Vcldman1  '  and  Moses^  *  .  All  of 

these  schemes  consist  of  coupling  a  system  of 
elliptic  equations  for  the  Inviscid  flow  to 
parabolic  equations  for  the  boundary  layer. 
However,  all  of  these  studies  neglect  the  fact 
that  the  parabolic  boundary  layer  equations  will 
not  give  accurate  resultc-for  a  thick  viscous 
layer  in  internal  flows ^  .  Higher  order 

boundary  layer  equations  are  treated  in  reference 
1  and  2  in  a  more  complete  manner  than  found  in 
most  of  the  literature  on  viscous- inviscid 
interaction.  However,  there  is  no  doubt  that  for 
"thick  boundary  layers,"  full  Navler-Stokes  or 
parabolized  Navier-Stokes  equations  are  far 
better  tools  for  treating  the  viscous  layer.  It 
is  our  goal  in  this  section  to  clarify  the 
viscous- inyiscid  interaction. 

For  thfc  purpose  of  this  study,  the  XLR-134 
nozzle,  which  we  computed  to  have  a  fairly  thick 
boundary  layer,  was  chosen.  At  each  axial 
station  flow  values  were  obtained  in  a  direction 
perpendicular  to  the  wall  (boundary  layer 
coordinate)  using  the  VNAP2  code.  Both  viscous 


fable  3:  Corperlton  0/  Ihrust  Calculations  for  1DK.  VNAP2(MS)  and  PKS  Codes  for  Viscous  flow  of  a  Perfect  Cit, 
Adiabatic  Wall 


lor/iiH 

VNAP2 

PNS 

NOZZLE 

HASS 

flow 

KAie, 

Ibm/sec 

THRUST, 

Ibf 

1 

CMOS 

AXIAL  x  CFL 

.  RADIAL _ 

l  SOftlC  TRANSONIC 
SUPERSONIC  SUPERSONIC 

HASS 

FLOW 

RATE 

lbm/*tc 

THRUST, 

Ibf 

ISP. 

aec.tbf 
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Table  4:  Subsonic  layer  Thickest  Predicted  by  VNAP2,  IlH  end  PNS  Codes 
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Engine 
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(BIN) 
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Subsonic  layer/ 
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loundiry  layer 

ft  (inches) 
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<S  (inches) 

0.995 
(BIN) 

Points  Across  the 
Subsonic  layer/ 

Points  across  the 
Boundary  layer 

6  (inches) 

Subsonic 

mu 

IlM 

PNS* 

VNAP2 

ILH 

PNS 

VNAP2* 

IlH 

13 

VNAP2 

UK 

PNS 

■ 

2.13 

3/14 

45/142 

5/16 

.0178 

.0072 

.155 

2.4J 

5/15 

jm 

4/13 

.0527 

.0513 

.2362 

2.02 

1/5 

58/126 

3/15 

.0071 

.034? 

.0551 

2.22 

4/18 

5/14 

.0560 

.0424 

.1541 

_ . 

*  Nurfcer  of  points  In  the  radial  direction  for  ASE  Is  95  and  for  "LR134  Is  01  with  more  points  dost  to  the  wall* 
This  is  an  approximation  for  VNAP2  and  PNS  Codes. 


and  invlscld  solutions  were  obtained  for 
comparison.  The  difference  between  the  axial  and 
total  velocities  for  the  inviscid  and  viscous 
solutions  are  depicted  in  Figure  1.  Also  shown 
on  Figure  1  is  the  displacement  thickness,  5*t 
as  computed  by  the  TDK/BLM  code.  In  classical 
TSL  theory,  the  effect  of  the  boundary  layer  on 
the  outer  or  core  flow  is  taken  into  account  by 
displacing  the  wall  in  the  normal  direction  by 
6*.  The  outer  flow  is  then  recomputed  using 
the  displaced  wall  and  thus  reflects  the 
interaction  of  the  TSL  on  the  core  flow.  The 
different  between  these  two  inviscid 
caiculal  is  (from  TDK)  divided  by  tho  original 
Inviscid  velocity  are  shown  on  Figure  1  and 
labeled  as  the  TSL  asymptote.  As  is  clearly 
evident  in  Figure  1,  the  fractional  velocity 
differences  approach  a  higher  asymptote  than 
predicted  by  TSL  theory  and  at  a  much  farther 
distance  from  the  wall. 


The  S3me  results  for  ASE  nozzle  are  shown  in 
Figure  2.  From  the  boundary  layer  code  6  - 

1.5  Inches  at  X  -  60  inches  (this  is  about  /I  of 
the  flowfleld,  i.e.,  6/ r  -  .07).  The  boundary 
layer  thickness  relative  to  the  local  transverse 
curvature  (S/t)  is  about  half  of  the  XLR134 
value.  So  it  is  anticipated  that  the  boundary 
layer  equations  yield  more  accurate  results  for 
the  ASE  nozzle.  ?rom  Figure  2,  it  can  be  seen 
that  at  6  -  1.5  inches  the  value  of 

(Ulnv'Uvl*)/Ulnv  ls  chan8ln8  very  slowly  «* 
opposed  to  the  region  closer  to  the  wall,  meaning 

that  the  viscous  efects  are  almost  dissipated  at 
the  1.5  inches  outward  from  the  wall.  Howevrr, 
there  is  still  some  viscous  effects  or 
interaction  for  distances  larger  than  1.5  inches. 

The  VNAP2  and  TDK  code  solutions  for  inviscid 
flow  of  a  perfect  gas  yield  the  edge  condition 
for  solving  the  boundary  laysr  aquations. 


Boundary  layer  equations  were  solved  for  oil 
cases  using  the  two  sets  of  edge  conditions  and 
the  results  are  indicated  In  Table  5.  This  table 
reveals  that  for  the  boundary  layer  losses,  the 
agreement  between  the  two  calculations  are  not 
good.  Since  the  boundary  layer  code  used  for 
both  calculations  was  the  same,  the  differences 
are  attributed  to  the  Invlscld  edge  conditions. 
Both  the  mesh  resolution  and  accuracy  of  the 
Method  of  Characteristic  calculation  are 
considered  superior  to  those  of  the  VNAP2  code. 

To  examine  the  magnitude  of  the  streamwlse 
diffusion  terns  In  Navler-Stokes  equations,  the 
full  NS  equations  are  compared  against  NS 
equations  without  the  streamwlse  diffusion  terms 
for  one  engine,  the  ASE.  Table  6  compares  the 
thrust  and  nass  flow  rate  for  both  cases,  ond 
Figure  3  Indicates  the  Mach  number  and  pressure 
at  the  exit  plane  for  both  cases.  Excellent 
agreement  Is  obse  ved.  Figure  3  Indicates  n 
severe  pressure  giadlent  at  the  wall  throughout 
the  boundary  layer  and  mild  pressure  gradient 
for  the  invlscld  portion  of  the  flow.  This  is 
compatible  with  the  results  in  Reference  1  and  2 
where  the  second  order  boundary  layer  solution 
reveals  this  pressure  gradient  across  the  shear 
layer.  These  sample  results  Indicate  that  the 
magnitude  of  streamwlse  diffusion  terms  In 
Navler-Stokes  equations  for  supersonic  flows, 
compared  to  other  terms,  are  very  small,  and  that 
they  can  be  eliminated  from  the  NS  equations. 

The  performance  results  with  axial  diffusion  terra 
is  closer  to  the  TDK/BLM  results  (see  Table  3). 
However,  the  difference  between  the  1  's  In 
Table  6  can  be  due  to  the  numerical  scRerae  rather 


than  the  axial  diffusion  terras.  This  has  been  a 
conclusion  of  studies  presently  being  performed 
at  Software  and  Engineering  Associate  Inc. 


Table  5:  Comparison  of  tha  Boundary  Layer  Losses 
Between  TDK/BLM  and  VNAP2/BLH  for 
Adiabatic  Wall 
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X#p,  lbf.  sec./lbm. 

TOK/BLH 

VNAP2/BLH 

\  Difference 
(TDK-VNAP2)/TDK 

SSME 

2.057 

2.184 

-  6.17 

Hughes  300:1 

3.43 

2.696 

15.37 

ASE 

4.789 

*  5.614 

•17.227 

XLR-134 

5.635 

6.854 

•21.633 

Tabla  4:  Hast  Flow  Kata,  Thrust  and  I  for  tha  ASE 
In  th*  fratanca  and  Abatnct  of  Axial 
Diffusion  Taras 
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-p  lb. 

*P  lba 

49.5* 

22,121 

446.51 

49.4) 

22,101 

445.15  ! 
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«» *on  betweon  the  core  flow  and  vincoun  layer  predicted 
ruil  NS  and  Th}.n  Shear  Layer  approximat Ion  por  XJ/R-13J1  llozzle. 
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Figure. <?.  Interaction  between  the  core  ''low  and  the  viscous  layer  predicted 
by  the  full  US  and  Thin  Shear  T,ayer  approximation  for  ASF.  Nozzle. 


COMPARISON  OF  PRESSURE  A T  THE  EXIT  PLANE 
FOR  VISCOUS  FLOW  OF  A  PERFECT  OAS 
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COMPARISON  OF  MACH  NUMBER  AT  THE  EXIT  PLANE 
FOR  VISCOUS  FLOW  OF  fl  PERFECT  CAS 


o  ywz  Mira  ixi*  iimwUN  o”Vi**  uimour  in*  ii/fioioi 


0.00  0.20"  0.40  O.w  0.00  |.U  t'SOl  «•«  «•»  ^ 

rwuwc.  run  I  I  5S2 ***** —  - 

PS  jure  3.  Conpirlsons  of  exit  plane  Mach  number  and  pressure  pronies  between  KG  solutions 
I’.h  ’ml  without  ntrcamwlsc  diffusion  For  ACS  Nozzle. 


vi.  awrutsTON 


The  effects  of  various  methods  of  computing 
the  wall  shear  layer  In  propulsive  nozzle*  has 
been  studied  for  relatively  thick  boundary 
layers.  It  was  found  that  for  nozzle  flows  where 
the  computed  boundary  layer  thickness  exceeds  10% 
of  the  local  nozzle  radius  that  calculations 
based  on  the  traditional  Prandtl  thin  boundary 
layer  assumptions  are  not  of  sufficient  accuracy 
for  detailed  propulsion  studies.  No  firm 
criteria  or  methods  have  been  developed  to 
estimate  the  magnitude  of  the  error  associated 
with  the  use  of  the  thin  shear  layer  assumptions 
The  work  that  this  report  Is  based  on  is  still  in 
progress  and  It  Is  hoped  that  more  definitive 
results  will  be  generated. 
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ABSTRACT 

A  Parabolized  Navier-Stokes  Code, 
VISPER  (Viscous  Performance  Code  for  Nozzle 
Flow  with  Finite  Rate  Chemistry)  has  been 
developed  which  calculates  the  internal  flow  of 
turbulent  and  non-equilibrium  reacting  gases. 
The  technique  of  Richardson  Extrapolation  is 
applied  to  the  1st  order  Beam-Warming 
scheme  to  provide  (i)  estimates  for  the  local 
error  at  each  marching  step,  (ii)  automatic  step 
size  variation,  and,  (iii)  extension  to  second 
order.  The  resulting  numerical  procedure  is 
more  stable  than  the  second  order 
Beam-Warming  method  and  has  the  efficiency 
of  an  atuoma;::  step  size  control  method.  Two 
turbulence  models,  mixing  length  and  k~e,  are 
used  to  resolve  the  wall  shear  layer.  The  results 
from  this  code  are  compared  against  existing 
experimental  data  for  supersonic  combustion 
and  rocket  nozzle  flows.  For  the  latter, 
comparisons  are  also  made  to  predictions  from 
classical  inviscid/boundary  layer  methods.  A 
measure  of  the  interaction  between  the  core 
flow/boundary  layer  is  obtained  and,  at  the 
same  time,  the  extent  of  the  validity  of  the 
classical  method  of  calculation  is  revealed. 


NOMENCLATURE 

c  mass  fraction 

Cp  specific  heat 

h  species  enthalpy 

H  total  enthalpy  (h  +  V2/2) 

J  Jacobian 

p  pressure 

Pf  prandtl  number 

q  heat  flux 

R  gas  constant 

r  transverse  curvature 
T  temperature 

V  velocity 

v  axial  and  normal  components  of  velocity 

x  axial  direction 

y  normal  direction 


Creek  Symbols 

£,  r\  transformed  coordinates  for  x  and  y 
co.  species  production  term 

0)  Vigneron  pressure  splitting  term 

(X  viscosity 

p  density 

t  shear  stress 


Superscripts 

*  reference  value 

-  contravariant 

i  index  in  £  direction  (axial) 

-  dimensional  quantity 
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Subscripts 

i  species 

j  index  in  T)  direction  (radial) 

L  laminar 

r  reference  value 

T,t  turbulent 

x,  £  in  axial  or  transformed  axial  direction 

y,  r\  in  normal  or  transformed  normal 

direcaon 
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1 


introduction 

In  recent  years,  there  has  been  a 
substantial  amount  of  research  in  computing 
steady,  viscous,  internal  flow,  due  to  its 
importance  in  the  design  of  components  such  as 
gas  turbine  and  ramjet  combustors, 
turbo-machinery,  inlet  ducts,  and  rocket  engine 
nozzles.  The  equations  that  govern  such  flows 
arc  the  Navicr-Stokc s  equations.  Under  certain 
conditions  these  flowfields  can  be  accurately 
predicted  with  a  simplified  version  of  the  full 
Navicr-Stokcs  equations.  The  Parabolized 
Navicr-Stokcs  equations  represent  one  such 
class  of  simplifications. 

Parabolized  Navier-Stokes  (PNS) 
equations  arc  a  subset  of  the  Navier-Stokes 
(NS)  equations  which  arc  valid  for 
predominently  supersonic  flow  with  subsonic 
shear  layers,  provided  no  streamwise  separation 
occurs1/.  The  PNS  equations  neglect  the 
streamwise  diffusion  term  which,  along  with 
special  treatment  of  the  subsonic  region  of  the 
boundary  layer,  removes  the  spatial  ellipticity 
from  the  steady  fonn  of  the  equations  and 
permits  a  solution  using  a  streamwise.  marching 
computational  technique.  Although  the  PNS 
models  were  developed  in  the  eorfy  1960's,  they 
were  not  widely  used  until  the  1970's  and 
1980's4 -9. 

As  the  interest  in  understanding  of  high 
speed  flows  and  performance  of  propulsive 
nozzles  increases,  there  is  need  for  development 
of  more  precise  computational  techniques.  It 
has  been  shown  that  the  adequacy  of  the 
standard  JANNAF  (Joint  Army,  Navy,  NASA, 
Air  Force)  liquid  rocket  performance  prediction 
methodology10  is  questionable  for  high 
expansion  ratio  nozzles  and/or  nozzles  with 
thick  viscous  layer.  The  JANNAF  procedure 
uncouples  the  inviscid  flow  and  viscous  layer 
and  loses  validity  when  there  are  extensive 
interactions  between  the  viscous  layer  and  the 
core  flow. 

The  PNS  equations  arc  integrated 
throughout  the  viscous  and  inviscid  regions  of 
the  flow,  This  procedure  eliminates  the  need  tc 
specify  the  edge  conditions  in  matching 
boundary  layer  and  inviscid  solutions,  i.e.  the 
conventional  inviscid-viscous  interaction. 
Popular  algorithms  for  solving  PNS  equations 
arc  those  by  Briley  and  McDonald11  and  Beam 
and  Warming12,  and  more  recently,  upwind 
algorithms41*4^. 


applied  extensively  to  boundary  layer  equations 
by  Keller  and  Ccbcci  36,  BlottncN7,  and 
others39,  cither  to  raise  the  order  of  the 
numerical  scheme  or  to  do  accuracy  study. 
More  recently,  Richardson  Extrapolation  has 
been  applied  to  Multigrid  Methods3*,40,  to 
provide  local  error  estimates  and  extension  to 
higher  order.  The  use  of  Richardson 
Extrapolation  here  differs  from  that  of 
references  36,37  in  two  ways.  Firstly,  the 
technique  is  applied  to  the  marching  direction 
and  not  the  radial  direction,  to  obtain  the  local 
error  and  extension  to  second  order.  Secondly, 
the  information  from  the  local  error  is  utilized 
to  provide  automatic  step  size  control  in  the 
marching  direction. 

The  resulting  algorithm  has  been 
implemented  in  a  production  code  for  liquid 
rocket  performance,  called  VISPER,  developed 
for  AFAL.  The  predictions  from  VISPER  are 
compared  to  experimental  data  as  well  as  to  the 
established  JANNAF  code,  TDK/BLM/ 
MABL29.  For  nozzles  known  to  have  a  thick 
viscous  layer,  comparing  VISPER  results  to 
classical  inviscid/boundary  layer  calculations 
yields  a  measure  of  the  interaction  between  the 
core  flow  and  the  boundary  layer. 
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For  body  fitted  coordinates  (5,q),  with  5 
«  £(x,r)  and  x\  *  q(x,r),  the  PNS  equations  for 
axisymmetric  flow  in  conservation  form  are: 


d 

3? 


,SXE  +  V 


+  3n 


|f  r.;'B+Ev)  +  71r(F+Fv)j 
+  ?  =  0  (1) 


where 


E  *  (pu,  pu^p,  puv,  pHu) 

EV 

F  «  (pv,  puv,  pu*+p,  pHv) 

FV  "C'VVWP 


Richardson  Extrapolation  is  a  generic 
numerical  procedure  that  can  be  applied  to  any 
computational  scheme  for  which  there  is  an 
asymptotic  expansion  of  the  local  truncation 
error  as  the  mesh  is  infinitely  refined27.  This 
extrapolation  technique  can  be  used  to  provide 
estimates  for  local  error  as  well  as  a  way  to 
extend  a  numerical  scheme  to  higher  order. 
Richardson  Extrapolation  has  been  used 
successfully  in  the  numerical  solution  of 
ordinary  differential  equation,  and  numerical 

Quadrature,  to  name  a  few  27.  For  fluid 
ynamic  applications  this  technique  has  been 


G  «  (0,  0,  -p, 

*x  Sr 

J  « 

Txx  s-^3Vx-3Vr) 

V  °-'l(Vr4Vx) 
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V  — l‘(3V--ty'x) 

qx  *  -  Pr  \\ 
qr  tyV 

For  planar  flow,  simply  remove  r  and  G 
from  equation  (1).  In  the  above  equations,  the 
viscosity,  ji,  represents  the  effective  viscosity, 
i.c.,  for  laminar  flow,  \i  is  the  laminar  viscosity, 
while  for  turbulent  flow,  p  is  the  sum  of  the 
laminar  and  turbulent  viscosity.  Similarly,  for 
turbulent  flow,  the  effective  thermal 
conductivity  is  represented  by 


The  thermodynamic  properties  are 
related  to  the  species  mass  fractions,  c.,  as 

follows: 

h  -iCih.cn 


rg 

the  summation  is  taken  over  the  number  of 
species. 

The  species  continuity  equations  for 
axisymmctric  flow  in  transformed  coordinates 
are: 


The  following  non-dimensionalization  is 
implemented 

x  »  xjt* 

h  «  RAi2 

r  «  r/x* 

P  -  p/pruj 

U  «  U/Uf 

\x  «  p/r*prur 

v  -  v/uf 

to.  m  coirVprur  , 

T  a  T/Tr 

and  consequently  the  form  of  equations  (1)  and 
(2)  remains  unchanged. 

FINITE  DIFFERENCE-EQUATIONS 

The  numerical  scheme  used  for  equation 
(1)  is  basically  the  linearized  block  implicit 
method  of  Briley-McDonald**  or 
Beam-Warming*2  The  present  implementation 
of  the  finite  difference  scheme  differs  from 
previous  PNS  solvers  in  two  major  ways.  First, 
the  vector  (p,  u,  v,  h)  is  used  as  dependent 
variables,  instead  of  the  more  common  vector 
(p,  pu,  pv,  pH).  Second,  the  diffusion  terms  are 
treated  using  second  derivatives  explicitly.  The 
advantages  of  the  present  approach  are  that  the 
tridiagonal  block  matrices  do  not  need 
conditioning  and  no  artificial  damping,  either 
implicit  or  explicit,  is  necessary.  The  geometric 
conservation  law15  is  observed  in  differencing 
the  fluxes. 


~3ci  -3ci  a 

Pu3r+pv35T3n 


+ 


Rr 


+  COj  + 


(2) 


The  species  continuity  equations, 
equation  (2),  are  finite-difference  using  the 
implicit  Euler  method  at  the  forward  station, 
once  a  solution  for  the  fluid  equations  is 
obtained  there.  This  decoupled  approach 
between  the  fluid  and  the  chemistry  cuts  down 
on  the  bandwidth  of  the  block  tridiagonal 
matrices,  while  does  not  seem  to  affect  the 
accuracy  of  the  scheme. 


Unity  Lewis  number  is  assumed  in 
equation  (2)  to  simplify  the  coding  and  should 
not  significantly  affect  the  accuracy  of  the 
results  for  the  problems  considered  here. 

For  planar  flow,  simply  remove  the  j 

term  in  equation  (2).  In  equation  (2),  u  and  v 
denote  the  contravariant  velocities, 


The  Vigneron  sublayer  model13,*4  is 
used  to  suppress  the  ellipticity  of  the  subsonic 
region  of  the  boundary  layer.  More  details  on 
the  numerics  are  available  in  reference  35,  as 
this  paper  is  oriented  toward  applications. 

Equation  (1)  is  rewritten  is: 


u  “^x  +  v^r  +G-0  (3) 

V  -  UT1X  +  VTlr 

b)  -  species  production  terms. 


3 


where 


.rV  +  V_fr 

•r - J - p 

■x<0,  ^x(l— co)P>  ^r(l-co)p.  0) 

s  parameter  in  Vigneron 
Sublayer  model, 


l  1  +  OH)M^  J 

0.7  £  o  £  0.85 

Mj:  *  Mach  Number  in  the  £ 

direction 

tix(E+Ev)  +  Tir(F+Fv) 

F  =  r - t - 


In  the  Beam-Warming  scheme,  the  flux 

vectors  E,  F,  G  need  to  be  linearized  about  the 
present  station.  Each  of  the  flux  vectors  can  be 
represented  as  of  the  sum  of  terms  of  the  form 
fW,  where  f  is  a  grid  related  quantity  such  as 
S  r/J.  r|x  r/J,  or  1/J,  and  W  is  a  flux  vector, 

such  as  E,  F,  or  G.  Denoting  by  j  the  index  for 
r|t  and  i  the  index  for  5,  then, 


Differences  in  the  x\  direction  are  eva'>  ated  as 
central  differences  at  the  forward  point.  When 
linearization  Of  the  form  (4)  asi  (5)  are  applied 
to  equation  ^3),  a  tridiagonm  block  structure 
results  which  can  be  solved  efficiently  for  the 
increment, 

'*0  »  Oi  +  1  — 

Vj  Vj. 

Once  a  solution  for  the  fluid  dynamic 
equations  (1)  is  obtained  at  the  forward  station, 
the  species  continuity  equations  can  be  solved 
at  the  same  station,  using  the  implicit  Euler 
scheme. 


At  the  solid  wall,  the  following 
bour  dary  conditions  arc  implemented: 

u  =  0  v  =  0 

3p/3n  «  0,  (3h/9r.  *  0  or  T  *  Tw). 

At  the  centerline,  symmetric  conditions  arc 
implemented: 

9u  _  rt  _  dp  _  dh  _  n 


where  n  is  the  norma:  vector  to  the  wall  and/or 
axis. 


-(fW) 


=Af)j  w!+1+fj+1AW! 

Using  Taylor’s  expansion  to  first  order. 
AWj  =W!+1-W! 

3W  An 

-31J4<5j 

(5> 

+  3^3ri(AQ)j 


(p,  U,  v,  h) 


Mixing  length16,17  and  k~e  models11,19 
for  the  effective  viscosity  have  been 
incorporated  into  the  VISPER  code.  An 
extensive  discussion  *bout  the  k~e  model  can 
be  found  in  Ref.  20.  The  method  used  to  obtain 
the  starting  conditions  for  the  model  were 
based  on  the  philosophy  that  the  starting 
solution  must  at  least  satisfy  the  approximate 
governing  equations  for  k  and  e.  To  obtain 
profiles  for  k  and  e  a  local  similarity  solution  of 
the  mean-flow  governing  equations  was 
obtained  using  the  mixing  length  model  at  the 
starting  marching  step21. 

This  similarity  solution  included 
approximate  profiles  for  the  eddy  viscosity  and 
the  ReynoL  s  shear  stresses.  The  initial  profiles 
for  k  and  e  could  then  be  evaluated. 

Some  uncertainty  regarding  the  nature  of 
turbulent  flow  for  low  values  of  the  local 
Reynolds  number  still  exists.  For  example,  the 
no-slip  boundary  condition  for  flow  over  a 
solid  surface  ensures  that  viscous  effects  will  be 
dominant  in  the  immediate  vicinity  of  the 
surface,  leading  to  low  local  Reynolds  numbers. 
Jones  &  Launder19  extended  the  k~e  model  to 
the  wall  region  without  using  a  near  wall 
analysis.  However,  their  model  has  not  yet 
been  thoroughly  verified. 
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Patankar  and  Spalding22  have  applied 
the  k~e  model,  incorporating  a  Couette  flow 
analysis  for  the  near  wall  region  and  using  a 
modified  van  Driest  formula23.  Arora,  ct  il2*,23 
have  compared  boundary  layer  solutions 
obtained  with  near  wall  models  to  several  sets 
of  experimental  data,  including  planar  and 
axisymmctric,  incompressible  and  compressible, 
laminar  and  turbulent,  reacting  and 
non-reacting  and  both  subsonic  and  supersonic 
flows.  Their  results  indicate  that  the  modified 
van  Driest  formula  suggested  by  Ccbcci  and 
Chang24  is  suitable  for  a  wide  range  of 
boundary  layer  flows.  Thus,  this  methodology 
was  adopted  for  the  present  effort. 

ARTIFICIAL  DISSIPATION 

In  all  the  test  cases  presented  in  Ms 
paper,  no  artificial  damping  is  used.  This  added 
stability  is  attributed  to  the  way  the  numerics 
are  implemented,  as  explained  in  the  section  on 
Finite  Difference  Equations.  Although  the 
usual  4th  order  explicit  and  2nd  order  implicit 
smoothing  terms  for  PNS  solvers1  have  been 
included  into  the  code,  it  has  been  found  that 
activating  either  of  them  is  detrimental  to  the 
accuracy  of  the  code.  This  fact  is  especially 
noteworthy  since  the  4th  order  smoothing  term 
does  not  affect  the  formal  accuracy  of  a  lower 
order  method. 

STEP  CONTROL  AND  RICHARDSON 

EXTRAPOLATION  VERIFIED 

The  technique  of  Richardson 
extrapolation27  has  been  applied  to  provide  an 
elaborate  step  size  control,  and  at  the  same 
time,  extend  the  numerical  scheme  to  second 
order  in  the  marching  direction.  Briefly 
speaking,  the  procedure  works  as  follows.  At 
each  marching  step  from  £  to  $  +  A£,  the 
solutions  at  £  +  sre  obtained  two  ways.  The 
first  way  takes  a  step  of  size  A£  and,  the  second 
takes  2  steps  of  size  A^/2.  The  results  at  £  + 
obtained  with  one  step  a;id  2  steps  are  then 
compared  to  obtain  the  local  error.  If  the  local 
error  exceeds  a  preset  limit,  then  the  step  is 
halved  and  the  procedure  repeated.  If  the  local 
error  is  within  the  accepted  limit,  then  the 
solution  at  ^  +  A£  is  extrapolated  to  second 
order  in  the  marching  direction  as 


where  f  denotes  any  dependent  variables,  and 
f\  f**  denote  ihe  solution  at  ^  +  A£  obtained  in 
one,  2  steps,  respectively. 

The  procedure  works  quite  well  for  this 
PNS  application.  Figure  1  shows  a  typical  step 
size  variation  as  a  function  of  axial  distance  for 
nozzle  flows.  In  this  flow  situation,  the  region 
of  highest  gradient  is  near  the  throat,  which 
requires  a  small  step  size,  and  as  the  flow 
expands  through  ihe  roiule,  a  larger  step  size 
can  be  used. 


om 

A  simple  algebraic  grid,  based  on 
transfinitc  interpolation28,  is  used  in  all  the  test 
cases  presented  in  this  paper.  A  sample  grid  is 
shown  in  Figure  2. 


FIN1TE.RATE.CHEMISTRY 

The  basic  finite  rate  chemical  kinetics 
capability  was  taken  from  the  TDK29  code. 
This  includes  the  reaction  rate  processor  and  the 
species  net  production  rate  evaluator. 
Capability  for  both  Arrhenius  and 
Landau*  'Viler  rate  data  forms  are  allowed. 
Both  second  and  third  order  reactions  are 
ireated  for  the  latter,  specific  3rd  body  reaction 
rate  ratios  can  be  input  to  the  code. 

TRANSONICLSOLUTION 

The  initial  start  line  for  the  PNS  is 
obtained  by  taking  the  inviscid  transonic 
analysis  from  the  Two-Dimensional  Kinetics 
Code  (TDK)29  and  attaching  a  boundary  fo^er 
next  to  the  wall.  The  above  transonic  solution 
is  approximate  and  based  on  the  perturbation 
method.  This  approach  is  justifiable  in  view  of 
the  very  thin  boundaiy  layer  at  the  throat.  One 
other  advantage  in  using  a  TDK  based  transonic 
solution  is  that  a  direct  comparison  with 
TDK/BLM29  (TDK/Boundary  Layer  Module) 
and/or  TDK/MABL  (TDK/Mass  Addition 
Boundary  Layer)  codes  can  be  made. 

RESULTSj^NDJBISCUS^IQN 

Three  test  cases  were  selected  to  verify 
the  code.  The  first  test  case  is  for  verification 
of  the  combustion  and  turbulence  models  and 
the  other  two  test  cases  deal  with  the  reacting 
flow  in  nozzles  where  the  viscous  layers  are 
known  to  be  thick. 

1.  Burrows  &  Kurlov30  supersonic 
combustion  to  verify  the 
combustion  and  turbulence 
models  in  the  code. 

2.  SBS-1 A  spacecraft  engine. 

3.  XLR-I34  engine. 

The  specifications  of  the  two  liquid  rocket 
engines  are  given  in  Tables  1  and  2,  together 
with  the  boundary  layer  thickness  computed 
using  a  classical  thin  boundary  layer  approach. 
It  can  be  seen  that  the  viscous  layer  is  fairly 
thick  in  these  nozzles  as  predicted  by  the 
boundary  layer  equations  and  for  these  cases  the 
extent  of  the  interaction  is  rather  significant10. 
In  these  cases  the  equations  should  be 
integrated  through  the  entire  flowficld  rather 
than  conventional  inviscid-viscous  interaction. 
So,  these  two  test  cases  seem  to  be  suitable  for 
these  studies. 
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Test  Case  1: 


Supersonic  Combustion^ 

This  test  case  was  employed  to  verify 
the  validity  of  the  combustion  and  turbulence 
models  within  the  PNS  code.  Figure  3  is  the 
schematic  of  the  test  ease  and  Figure  4 
cc-mpares  the  experimental  data  for  the  mass 
fraction  of  H20  with  what  is  predicted  by  the 
PNS  code  using  mixing  length  and  k*-£ 
turbulence  models.  From  these  results  it  can  be 
concluded  that  accurate  predictions  of  species 
mass  fractions  are  obtained  with  VISPER,  and 
that  little  difference  is  seen  between  the  two 
turbulence  models. 


Test  Case  2:  SBSr.lA.Spflv.ccrafi Engine 


A  Hughes  Aircraft  Corporation  (MAC) 
small  thruster,  IIE-54,  was  used  to  validate  the 
theoretical  performance  prediction  of  the  PNS 
code.  This  engine  was  selected  for  our  study 
for  the  following  reasons.  First,  experimental 
data  are  readily  available,  as  the  engine  had 
been  studied  extensively  by  Kushida  ct  al  31. 
The  second  reason  is  the  flowficld  through  the 
nozzle  is  completely  laminar  resulting  in  a  very 
thick  boundary  layer.  The  previous  study31  was 
done  using  the  classical  inviscid-boundary 
layer  approach.  The  present  study,  using  the 
VISPER  code,  should  be  able  to  address  the 
adequacy  of  the  latter  approach. 


This  engine  is  a  small  station  keeper 
engine  using  decomposed  hydrazine  as 
propellant.  Hence,  all  of  the  PNS  calculations 
were  performed  using  3  species,  NH3,  N2,  and 
H2,  and  frozen  chemistry.  The  chamber 
operating  conditions  were  as  described  in 
Reference  31,  and  arc  repeated  here  as  Table  3. 
Six  chamber  pressures  ranging  from  25  psi  to 
132  psi  were  considered.  Performance 
predictions,  in  terms  of  the  thrust  coefficient, 
were  obtained  for  all  six  eases.  However,  the 
flowficld  results  arc  only  discussed  for  the  two 

extreme  cases,  namely  the  P„  =  25  psi,  and  P 
c  c 


For  the  first  ease,  with  Pc  =  25, 


=  132  psi. 

Figures  5  and  6  show  the  pressure  and  Mach 
number  contour  line.  The  thickness  of  the 
boundary  layer,  especially  the  sonic  line,  (M=l) 
can  be  seen  from  Figure  6.  Both  Figures  5  and 
6  show  very  clearly  a  compression  wave 
coming  Fom  the  attachment  point  in  the  nozzle 
contour  and  reflecting  off  the  centerline.  The 
magnitude  of  the  compression  wave  can  be  seen 
from  Figure  7,  which  shows  that  the  pressure 
along  the  centerline  increased  significantly.  It 
is  interesting  to  note  that  this  nozzle  was 
designed  as  a  boundary  layer  compensated  Rao 
nozzle.  The  fact  that  a  shock  wave  still  exists 
in  the  nozzle  indicates  that  this  classical  design 
method  may  not  be  applicable  to  flow  situations 
having  such  a  thick  boundary  layer.  Figures  8 


contours  for  the 


*  psv^urc  and 
Pc  =  132  ease. 


nmW 


The  boundary 


layer  is  still  very  thick,  although  not  as  thick  as 
for  the  Pc  «  25  case. 
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For  the  performance  prediction,  the 
thrust  coefficient,  Cp  defined  as  total  thrust 

divided  by  the  product  of  the  chamber  pressure 
and  the  throat  area,  were  obtained  from  the  6 
cases  using  the  PNS  code  and  then  compared  to 
the  results  of  Kushida3*  ct  al.  In  Reference  31, 
the  experimental  data  were  compared  against 
the  theoretical  prediction  of  the  JANNAF 
methodology  using  TDK/BLM29.  In  Figure  10 
the  results  of  Kushida  are  compared  against  the 
present  theoretical  performance  predictions.  It 
can  be  seen  in  Figure  10  that  the  TDK/BLM 
results  compare  surprisingly  well  with  the 
experimental  data,  considering  how  thick  the 
boundary  layers  are,  especially  in  the  low 
chamber  pressure  range.  However,  the  PNS 
results  are  in  better  agreement  with  the  data. 
Both  codes,  TDK/BLM  and  PNS,  show 
excellent  agreement  at  the  highest  chamber 
pressure,  which  is  understandable  given  the 
thinner  boundary  layer. 

Test  Case  3:  XLRt13-4,Nq«1c 

The  Aerojet  XLR-134  is  a  low  thrust 
cryogenic  engine  being  developed  for  the 
AFAL  (Air  Force  Astronautic  Laboratory).  It  is 
intended  to  provide  low  thrust  propulsion  for 
the  delivery  of  large  space  structures  to 
geosynchronous  orbit32,33-  The  nozzle  has  an 
exit  area  ratio  of  767.  The  injector  uses  a 
conventional  coaxial  element  in  which  the 
gaseous  fuel  flows  in  an  anulus  around  the 
liquid  oxidizer.  To  minimize  the  amount  of 
heat  transfer  from  the  engine  to  the  vehicle,  the 
entire  nozzle  is  regenerative  cooled.  It  should 
be  noted  that  the  laminar-turbulent  transition 
based  on  Rcq-  360  (Reynold’s  number  based  on 

the  momentum  thickness)  reveals  that  for  the 
given  operating  conditions,  the  flow  stays 
laminar  till  about  one  inch  before  the  throat. 
As  yet,  hot  firing  data  arc  not  available  for  this 
engine.  Thus,  strictly  speaking,  the  XLR-134 
is  not  a  validation  case.  However,  it  is  useful 
for  exploring  the  limits  of  the  boundary  layer 
method. 

The  calculations  for  this  engine  were 
performed  with  six  (6)  species  (H,  H2,  H20,  O, 
OH,  O2)  and  eight  (8)  reaction  rates.  Figure  1 1 
compares  the  pressure  along  the  wall  calculated 
by  the  VISPER  code,  against  the  TDK/MABL 
code.  In  this  nozzle  the  boundary  layer 
thickness,  (5.995/r)cxjt,  is  about  0.2,  which  can 

be  considered  a  fairly  thick  shear  layer. 

This  thickness  implies  that  the  viscous 

effects  can  extend  beyond  that  of  classical 

boundary  layer  predictions.  Furthermore,  the 

standard  methodology,  i.e.,  boundary 

layer/inviscid  flow  interaction,  can  undcrprcdict 

the  specific  impulse  I  (  =  thrust/mass  flow 
sp 

rate),  because,  in  general,  it  neglects  part  of  the 
viscous  effects  in  the  nozzle,  (see  Table  4). 
Figure  12  compares  the  wall  shear  between  the 
two  codes.  The  pressure  and  Mach  number 
contours  depicted  in  Figures  13  and  14  reveal 
that  the  shock  originating  from  the 
attachment  point  propagates  to  the  axis  and  then 
reflects. 
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Richardson  Extrapolation 
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Tabic  4  shows  the  run  time  and  I  for 

sp 
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first  half  of  the  table  is  for  the  results  with 
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supersonic  nozzle  flows  the  required  accuracy 
(<  1%)  is  difficult  to  attain.  It  has  been  shown 
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predictions  for  ISp10*34. 


CQIiCLLlDINQ-REMARKS 

The  Richardson  Extrapolation  technique 
has  been  successfully  implemented  for  a  PNS 
solver.  Its  advantages  arc  shown  both  in  the 
computer  run  time  and  accuracy.  This 
application  of  Richardson  Extrapolation  can  be 
extended  to  any  other  marching  flow  solver, 
such  as  time-dependent  Navier-Stokes  solver, 
the  only  restriction  is  the  flow  cannot  have 
strong  shock  wave. 

Th.  VISPER  code  has  shown  to  be  a 
viable  too!  for  studying  supersonic  nozzle 
flows,  especially  those  with  thick  viscous 
layers.  Excellent  agreement  with  experimental 
data  were  obtained  for  supersonic  combustion 
and  nozzle  flows.  With  respect  to  the 
interaction  of  core  flow/boundary  layer,  the 
classical  decoupled  approach  seems  to  be  valid 
to  nozzle  flows  with  a  Reynold's  number  (based 
on  throat  radius)  about  10,000  (test  case  2). 
However,  the  VISPER  code  provides  a  better 
alternative  to  the  decoupled  approach.  Not  only 
arc  the  results  more  accurate,  but  the  computer 
execution  time  is  much  less,  especially  if  a 
second  order  correction  is  performed  on  the 
core  fiow/boundary  layer  calculations  by 
adjusting  the  body  by  the  displacement 
thickness. 
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TNS  SOLUTION  OF  NON- EQUILIBRIUM 
REACT1NO  rLOW  IN  ROCKET  NOZZLES 

A.  L.  Dang*,  II.  Kehtarnavar**,  and  D.  E.  Coocs^ 

Software  end  Engineering  Associates,  Inc. 

Carson  City,  Nevada 

ABSTRACT 

A  Parabolized  Navler*5tokes  Code  (PNS)  which  calculates  the  Internal  flow  of  turbulent 
ami  tion-cqull Ibrlum  reacting  gasos  has  boon  dovolopod.  A  second  order  method  with  automatic 
step  size  control,  based  on  the  Beam ‘Warming  scheme  and  Richardson  extrapolation  technique, 
is  used  to  solve  the  flowflold  equations.  The  results  from  this  code  are  compared  against 
the  existing  JANNAF  performance  prediction  methodology  calculations.  The  rcsuRs  indicate 
very  good  agreement  and  reveals  the  advantages  of  the  Richardson  extropolat Ion  technique, 
both  in  terms  of  accuracy  and  computer  execution  time.  Two  turbulonco  models,  mixing  length 
and  k  -  <,  are  used  to  resolve  the  wall  shear  layer. 

INTRODUCTION 

Parabolized  Navier ‘Stokes  (PNS)  equations  jrej  a  subset  of  the  Navler-Stokes  (NS) 
equations  which  are  valid  for  supersonic  flows  ’  .  The  PNS  equations  neglect  the  strcnmwlsc 
diffusion  term  which,  along  with  special  treatment  of  the  subsonic  region  of  the  boundary 
layer,  removes  the  spatial  elliptlclty  from  the  steady  form  of  the  equations  and  permits  a 
solution  using  a  strennwlsc  marching  computational  technique.  Although  the  PNS  mode^SgWere 
developed  In  the  early  19G0's,  thoy  wore  not  widely  used  until  the  1970' s  and  1980's  *  . 

As  the  Interest  in  understanding  of  high  speed  flows  and  performance  of  propulsive 
nozzles  increases,  there  is  need  for  development  of  more  precise  compujgtlonal  techniques. 

It  has  been  shown  that  the  adequacy  of  the  standard  JANNAF  methodology  Is  questionable  for 
high  expansion  ratio  nozzlas.  The  current  procedure  consists  of  uncoupling  the  Invlsctd 
flow  and  viscous  layer  which  loses  validity  when  there  Is  extensive  interactions  between  the 
viscous  layer  and  the  core  flow. 

The  PNS  equations  are  integrated  throughout  the  viscous  and  inviscld  regions  of  the 
flow.  This  procedure  eliminates  the  need  to  specify  the  edge  conditions  in  matching 
boundary  layer  and  inviscld  solutions,  l.e.  the  conventional  invlswid. viscous  interaction. 
Furthermore,  the  PNS  equations  are  objjlned  by  Implicit  mothytjs.  Notable  algorithms  of  this 
type  arc  those  by  Briley  and  McDonald  *  and  Beam  and  Warming  , 

GOVERNING  EQUATIONS 

For  body  fitted  coordinates  (f,q),  with  (  -  $(x,r)  and  q  -  q(x,r),  the  PNS 
equations  for  ax l symmetric  flow  in  conservation  form  are: 
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For  planar  flow,  simply  remove  r  and  C  from  equation  (1).  In  the  above  equations*  the 
viscosity,  jj,  represents  the  effective  viscosity,  i.c.,  for  laminar  flow,  /j  is  the 
laminar  viscosity,  while  for  turbulent  flow,  /j  is  the  sura  of  the  laminar  and  turbulent 
viscosity.  Similarly,  for  turbulent  flow,  the  effective  thermal  conductivity  is  represented 
by 


Th«  thermodynamic  properties  ars  related  to  the  spscies  mass  fractions,  ,  as  follows: 

h  -  Ec  h(T)  cp  -  Ec  cp  (T)  R  *  Ec  R  (T) 

i  i  i  1  G  i  Cl 

The  species  continuity  equations  (unity  Lewis  number)  for  axisyrametric  flow  in  transformed 
cooidlnates  are: 


+  /V 


d  /j  9c^ 

7  (—  7-  n  )n 
dij  pr  3ij  xx 


J  v  3c ^  Ip 

+  --  ( —  ---  rj  )rj  +«+--- - t) 

in  pr  dn  r  r  i  r  pr  d?  r 
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For  planar  flow,  simply  remove  the  *  term  in  equation  (2).  In  aquation  (2), 


u 

-  uf  +  v( 
sx  r 

v  -  ur?x  +  vr?r 

-  species  production  terms 

The  following 

normalization  is 

implemented 
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and  consequently  equations  (1)  and  (2)  remain  exactly  the  same. 


FINITE  DIFFERENCE  EQUATIONS 
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The  numerical  scheme  used  for  equation  (1)  is  the  Beam-Warning  method,  while  ths 
scheme  used  for  eq^tjgn  (2)  Is  the  implicit  Euler.  Treatment  of  departure  solution  la  done 
using  the  VlgneronA  '  tcjjinlque.  The  flux  vectors  ore  linearized  according  to  the 
Geometric  Conservation  Law  . 


Thus,  equation  (1)  Is  rewritten  as: 


d  A  9  A  9  A  A 
?r  (E)  +  ;;<P)  +  7-  (F)  +  G  -  0 
d(  9( 

f  E  +  £_F  . 

where  E  -  r - j . P  P  -  r(0,  <x(l-w)p.  {r(l-u)p,  0) 

:  ”x<E+Ev>  +  VF+V  r.  G 

r  "  r . j .  °  ~  J 


(3) 


In  the  Bcara-Warralng  scheme,  the  flux  vectors  E,  F,  C  need  to  be 
present  station.  Each  of  the  flux  vectors  can  be  represented  as  of 
form  CU,  where  f  Is  a  grid  related  quantity  such  os  f  r/J ,  n  r/J  or 
vector,  such  as  E,  F,  or  C.  Denote  by  J  the  Index  for  rj*  and  1  the 


linearized  about  the 
the  sum  of  terms  of  ths 
1/J,  and  W  Is  a  flux 
index  for  Then, 
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A(fU)  -  (fW) 


(6f)j+1  t)j  +  fj+l  Auj 


Using  Taylor*#  expansion  to  fir#t  order: 


AWJ  "  WJ  1  #  Wj 


IS  6Qj +  IS;  aVj  -  hi  AQj +  IS; !;  <iQ)j 


Q  -  (p,  u.  V,  h) 


q  -!9  . 

fj  dn 


Vlicn  linearization  of  the  form  (A)  and  (5)  are  applied  to  equation  (3),  a  tridiagonal  block 
structure  results  which  can  be  solved  efficiently  for  the  increment, 

- «}“  •  i  • 

Once  a  solution  for  the  fluid  dynamic  equations  (1)  is  obtained  at  the  forward  station, 
tho  species  continuity  equations  can  be  solved  at  the  some  station,  using  the  implicit  Euler 
scheme . 

BOUNDARY  CONDITIONS 

At  the  solid  wall,  the  following  boundary  conditions  are  implemented: 


u  m  o  v  •  0 


ap/dn  -  o,  (ah/an  »  o  or  T  •  T  )• 

v 


At  the  centerline,  symmetric  conditions  are  implemented: 


au  dp  ah 

«  0  v  ■*  0  , 

an  an  an 


where  n  is  the  normal  vector  to  the  wall  and/or  axis. 


TURBULENT  MODELS 


Mixing  length  *  and  k-<  models  ’  for  the  effective  viscosity  have  been 
incorporated  into  the  PNS  code.  An  extensive  discussion  about  tho  k-e  model  can  be  found 
In  Ref.  (20).  The  method  used  to  obtain  the  starting  conditions  for  tho  k-<  model  were 
based  on  the  philosophy  that  the  starting  solution  must  at  least  satisfy  the  approximate 
governing  equations  for  k  and  i.  To  obtain  profiles  for  k  and  <  a  local  similarity 
solution  of  the  mean-flow  governing  equations  was  obtained  using  the  mixing  length  modol  at 
die  starting  marching  step 

This  solution  included  approximate  profiles  for  the  eddy  viscosity  and  the  Reynolds 
shear  stresses.  The  initial  profiles  for  k  and  t  could  then  be  evaluated. 

Some  uncertainty  regarding  the  nature  of  turbulence  flow  when  tho  local  Reynolds  number 
Is  low  still  exists.  For  example,  the  no-slip  boundary  condition  for  flow  over  a  solid 
surface  ensures  that  viscous  effects  will  be  dominant  in  the  immediate  vicinity  of  the 
surface,  leading  to  low  local  Reynolds  numbers.  Jones  A  Launder1  extended  the  k-<  model 
to  the  wall  region  without  using  a  near  wall  analysis.  Howovor,  their  model  has  not  yet 
beers  thoroughly  ver  i fled. 
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Tatankar  and  Spalding  have  applied  the  k-<  model,  incorporating  a  Cuuette  flow 
nn3iy5i*  ^or  C^e  near  wal*  ro6l°n  *nd  using  a  modified  van  Driest  formula^  .  Arora,  et 
nl  ’  have  compared  sovernl  sets  of  experimental  data,  including  planar  and  axlsymmetric , 
incompressible  and  compressible,  laminar  ond  turbulent,  reacting  and  non-reacting  and  both 
subsonic  and  supersonic  flows-,  Their  results  indicote  that  the  modified  van  Driest  formula 
*l*W&®*e*d  by  Cobeci  and  Chang  is  suitable  for  a  wide  range  of  boundary  layer  flows.  The 
came  methodology  has  been  used  in  this  work. 


I 


ARTIFICIAL  DISSIPATION 

In  all  the  test  cases  presented  In  this  paper,  no  artificial  damping  Is  used.  Although 
the  usual  6th  order  explicit  and  2nd  order  Implicit  smoothing  terms  for  PNS  solvers8  havs 
lu'cn  Included  into  the  code,  It  has  bcon  found  that  activating  either  of  them  Is  detrimental 
to  the  accuracy  of  the  code.  For  test  case  2  for  example,  using  Just  10\  of  the  recommended 
, (mount  for  t ho  explicit  smoothing  causes  the  I  to  rise  by  It,  which  Is  not  acceptable. 

Mils  fact  Is  especially  noteworthy  aince  the  order  smoothing  term  does  not  affect  the 
formal  accuracy  of  n  lower  order  method. 


STEP  CONTROL 

An  elaborate  step  control  procedure  has  bc(jtj  developed  for  the  PNS  solver.  This  method 
Is  based  on  Richardson  extrapolation  technique  which  has  been  used  successfully  In 
numerical  solutions  of  ordinary  differential  equations  and  other  numerical  calculations.  In 
order  to  make  use  of  the  Richardson  extrapolation  technique  for  the  present  application,  the 
PNS  Is  castcd  in  o  strictly  marching  procedure.  At  each  marching  step  from  £  to  (  f 
A*,  the  solutions  at  <  +  A(  are  obtained  two  ways.  The  first  way  takes  a  step  of  size 
A*  from  (  to  to  £  +  A(.  The  second  way  takes  2  steps  of  size  A(/2  from  {  to  ( 
i  a{/2  and  then  to  (  +  A{.  The  results  ac  (  +  A£  obtained  with  one  step  and  2 
steps  are  then  compared  to  determine  the  local  error,  which  in  turn  decides  on  the  next  step 
size.  The  procedure  works  quite  well  for  smooth  flow  situations,  but  requires  more  work  for 
flow  with  a  strong  shock. 


GRID 

A  simple  algebraic  grid,  based  on  transflnlte  interpolation^,  is  used  in  all  the  test 
cases  presented  In  this  paper.  A  sample  grid,  for  the  SSME,  is  shown  in  Figure  1.0. 

TRANSONIC  SOLUTION 

The  initial  start  line  for  the  IJtJS  is  obtained  by  taking  the  transonic  analysis  from  the 
Two-Dimensional  Kinetics  Code  (TDK)  and  attaching  a  boundary  layer  next  to  the  wall.  This 

approach  is  Justifiable  in  view  of  the  very  thin  boundary  layer  at  the  throat.  One  other 
advantage  In  using  a  TDK  based  transonic  solution  is  that  a  direct  comparison  with  TDK/BLM 
(TDK/Boundary  Layer  Module)  and/or  TDK/MABL  (TDK/Mass  Addition  Boundary  Layer)  codes  can  be 
made . 


RESULTS  AND  DISCUSSION 

Four  test  cases  were  selected  to  verify  the  code.  The  test  cases  arc; 
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1.  Burrows  &  Kurlov  supersonic  combustion  to  verify  the  combustion  and  nixing  model*  In 
the  code, 

2.  SBS-lA  spacecraft  engine. 

3.  XLR-136  engine. 

6.  SSME  (Space  Shuttle  Main  Engine), 

The  specifications  of  these  three  engines  are  given  in  Table  1, 

Test  Case  1;  Burrows  &  Kurlov  Supersonic  Combustion 

This  test  case  was  employed  to  verify  the  validity  of  the  turbulence  model*  within  tha 
PNS  code.  Figure  1.1  Is  the  schematic  of  the  test  case  and  Figure  1.2  compare*  the 
experimental  data  for  the  mass  fraction  of  H^O  with  what  is  predicted  by  the  PNS  code  u*ing 
raising  length  and  k-<  turbulence  models.  From  these  results  it  can  be  concluded  that  both 
turbulent  models  yield  fairly  satisfactory  results  for  supersonic  combustion. 

Test  Case  2:  SBS«1A  Spaceacraft  Engine 

A  Hughes  Aircraft  Corporation  (MAC)  small  thruster,  HE-56,  was  used  to  validate  the 
theoretical  performance  prediction  of  the  PNS  code.  This  engine  was  selected  for  our  study 
for  the  following  reasons.  First,  experiment‘d  data  are  readily  available,  as  the  engine 
had  been  .  ;udlod  extensively  by  Kushlda  et  alJ  .  The  second  reason  Is  the  flowfleld  through 
the  nozzle  Is  completely  laminar  from  the  throat  on,  resulting  ix\  a  very  thick  boundary 
layer.  The  previous  study3  wae  done  using  the  the  classical  invl*cid»boundery  layer 
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approach.  The  present  study,  using  the  PNS  cods,  should  bs  able  Co  sddrsss  tbs  sdsquscy  of 
the  latter  approach. 

This  engine  Is  a  small  scacion  keeper  engine  using  decomposed  hydrazine  as 
propellant.  Hence,  all  of  the  PNS  calculations  are  performed  using  3  species,  NH3,  N2,  and 
112.  and  frozen  chemistry.  The  chamber  operating  conditions  are  as  described  In  Reference 
31.  and  are  repented  here  ns  Table  2.  Six  chamber  pressures  ranging  from  25  psl  to  132  pal 
vc  ro  considered.  Performance  predictions.  In  terras  of  the  thrust  coefficient,  were  obtained 
for  nil  six  cases.  However,  the  flowficld  results  are  only  discussed  for  tht  two  extreme 
c.v:c*s,  nnraely  the  P  -  25  psl,  and  P  -  132  psl.  For  the  first  case,  with  P  -  25,  Figures 

2.1  and  2.2  show  the  pressure  and  Hach  number  contour  line.  The  thickness  of  the  boundary 
layer,  especially  the  sonic  line,  (M-l)  can  be  seen  from  Figure  2.2.  Both  Figures  2.1  and 

2.2  show  very  clearly  a  compression  wave  coming  from  tho  attachment  point  In  the  nozzle 
contour  and  reflecting  off  the  centerline.  The  magnitude  of  the  compression  wave  can  be 
seen  from  Figure  2.3,  which  shows  that  the  pressure  along  the  centerline  Increased 

signif icnntly.  It  Is  interesting  to  note  that  this  nozzle  was  designed  as  a  boundary  layer 
rompensnted  Rao  nozzle,  In  fact  by  SEA,  Inc.  for  Hughes  Aircraft  Corporation. 

The  fact  that  a  shock  wave  Is  still  in  the  nozzle  shows  that  this  classical  design 
raothod  rany  not  be  applicable  to  flow  situations  having  such  a  thick  boundary  layer.  Figures 
2.4  and  2.5  show  the  pressure  and  Hach  number  contours  for  the  P  -  132  case.  The  boundary 
layer  Is  still  very  thick,  although  not  as  much  as  the  Pc  ••  25  case. 

For  the  performance  prediction,  the  thrust  coefficient,  Cp,  defined  as  total  thrust 
divided  by  the  product  of  the  chamber  pressure  and  the  throat  area,  ^e  obtained  from  the  6 
cases  using  the  PNS  code  and  then  compared  to  the  results  of  Kushlda  et  al .  In  Reference 
31,  the  experimental  data~uere  compared  against  the  theoretical  prediction  of  the  JANNAF 
methodology  using  TDK/BLM  .  In  Figure  2.6  the  results  of  Kushlda  are  compared  against  the 
present  theoretical  performance  predictions.  It  can  be  seen  in  Figure  2.6  that  the  TDK/BLM 
results  compare  surprisingly  well  with  the  experimental  data,  considering  how  thick  the 
boundary  layers  are,  especially  in  the  low  chamber  pressure  range.  The  PNS  results, 
however,  are  In  excellent  agreement  with  the  data.  Both  codes,  TDK/BLM  and  PNS,  show 
excellent  agreement  at  the  highest  chamber  pressure,  which  correlates  well  with  the  thinner 
boundary  layer. 

Test  Case  3:  Xl.R-134  Nozzle 

The  Aerojet  XLR-134  Is  a  low  thrust  cryogenic  engine  being  developed  for  the  AFAL  (Air 
Force  Astronaut Ic  Laboratory).  It  is  intended  to  provide  Jgw^hrust  propulsion  for  the 
delivery  of  large  space  structures  to  geosynchronous  orbit  *  .  The  chamber  has  an  exit 

area  ratio  of  7.7.  The  Injector  uses  a  conventional  coaxial  element  In  which  the  gaseous 
fuel  flows  in  an  anulus  around  the  liquid  oxidizer.  The  injector  faceplate  Is  a  Jjjlck 
copper  section  to  provide  adequate  cooling  and  is  integral  with  the  injector  body  .  One  of 
the  major  considerations  in  nozzle  design  was  to  minimize  the  amount  of  heat  transfer  from 
the  engine  to  the  vehf  1 e.  For  this  reason,  regenerative  cooling  to  the  exit  of  the  nozzle 
vos  required.  It  should  be  noted  that  the  laminar* turbulent  transitions  based  on  Re* 
(Reynold's  number  based  on  the  momentum  thickness)  reveals  that  for  the  given  operating 
conditions,  the  flow  stays  laminar  till  about  one  inch  before  the  throat.  As  yet,  engine 
firing  data  is  not  available  for  this  engine. 

The  calculations  for  this  engine  are  performed  with  six  (6)  species  (H,  H«0,  0,  OH, 

0^)  and  eight  (8)  reaction  rotes.  Figure  3.1  compares  the  prossuro  along  the  wall 
calculated  by  the  PNS  couc,  against  the  TDK/MABL  code.  It  con  be  seen  that  PNS  predicts 
about  0.126%  more  loss  than  the  TDK/MABL  code.  In  this  nozzle  the  boundary  layer  thickness, 
($ . 995/r)^^ t ,  is  about  0.2,  which  can  be  considered  a  fairly  thick  shear  layer. 

Computer  Code  I  ,  lbf/ibm-sec 


TDK/HABL  457.45 

PNS  456.87 

This  thickness  implies  chat  the  viscous  effects  con  extend  beyond  that  of  classical 
boundary  layer  predictions.  Furthermore,  Che  standard  methodology,  *.c.,  boundary 
loycr/lnvlscld  flow  interaction,  can  overpredlcc  the  specific  impulse  because,  in  general, 
t liny  neglect  part  of  the  viscous  affects  in  the  nozzle.  Figure  3  2  compares  the  wall  shear 
heewoen  tha  two  codes.  Tha  pressure  and  Hach  number  plots  depleted  in  Figures  3.3  end  3.4 
reveals  that  tha  shock  originating  from  the  attachment  point  develops  to  the  axis  and  then 
reflects. 
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Test  Case  4: 


Space  Shuttle  Haln  Engine  (SSHE) 


The  Space  Shuttle  Haln  Engine  (SSHE),  which  Is  a  flight  engine,  utilizes  a  thrust 
chamber  that  Is  constructed  In  two  jgcijjons:  1)  a  Haln  Combustion  Chambor  <HCC)  assembly; 
anil  2)  an  expansion  nozzle  assembly  •  .  The  nozzle  assembly  Is  bolted  to  the  HCC  at  a 
nozzle  attach  flange.  The  expansion  ratio  at  the  attachment  position  Is  5:1.  The  nozzle 
expansion  assembly  extends  to  an  expansion  ratio  of  77.5:1. 

The  SSHE  uses  a  hydrogen  expander  cycle  with  oxygen  preheating  and  a  nozzle  that  Is 
logoncrativcly  cooled  with  hydrogen.  Nozzle  geometry  and  operating  conditions  for  the  SSHE 
current  design  nozzle  (109%  power  level)  were  obtained  from  NASA/MSFC. 

The  laminar- turbulent  transition  based  on  Re^  indicates  that  the  flow  In  the  chamber  and 
nozzle  Is  completely  turbulent.  The  boundary  layer  is  very  thLn  In  this  nozzle, 

(tf.995/r)  . c**0.0A ,  and  the  conventional  lnvlscld-vtscous  Interaction  yields  good  results 

for  the  performance  calculations.  The  PNS  calculations  for  the  SSHE  were  done  with  six  (6) 
species  (H,  ll^,  H^O,  0,  OH,  0^)  and  eight  (8)  reaction  rates. 


Figures  4.1  and  4.2  compare  the  pressure  and  total  wall  shear  along  tho  nozzle  wall 
between  the  TDK/HABL  and  PNS.  From  Figure  4.2  it  can  be  seen  that  PNS  predicts  more  los*  ts 
(higher  shear  at  the  wall). 


Computer  Codo 


I  lbf/lbnt.  see . 


TDK/HABL 

455, 

.96 

PNS 

451, 

.03 

The  dlf^grence  In  1^  's  ts  duo  to  the  fact  that  tho  chemistry  In  this  engine  is  in 
equilibrium  and  the  resulting  equations  become  stiff.  The  present  version  of  the  PNS  code 
does  not  properly  account  for  this  effect.  The  code  will  be  modified  to  handle  these  severe 
conditions.  Figures  4.3  or.  4.4  show  the  calculated  pressure  and  Hach  Number  contours.  The 
compression  wove  originating  from  the  attachment  point  can  be  clearly  seen. 

CONCLUDING  REMARKS 

The  results  obtained  In  this  paper,  although  still  preliminary,  are  very  encouraging. 

Hie  PNS  codo  has  shown  to  he  a  viable  tool  for  performance  prediction  over  a  wide  range  of 
conditions  for  a  variety  of  rocket  engines.  The  results  are  especially  good  for  flows  with 
thick  boundary  layers.  Some  problems  still  remain  to  be  investigated  and  resolved. 

For  example,  the  1%  underprediction  of  tho  SSHE,  and  for  nozzles  with  strong  shocks.  For 
che  latter,  an  upwind  scheme  may  bo  required. 
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figure  1.0.  Physical  Grids. 
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Proooure  Contour  for  Pc**132  for  tho  SRS-IA 
Spacecraft  Engine. 
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Figure  3.1,  Pressure  along  the  Figure  3*2. 
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Figure  3.4.  Maoh  Number  Contour#  for  XLR-131!  notslo 
turbulent  flow  with  chemistry. 
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Figure  4,3.  Pressure  Contours  for  SSME  noxzle 
turbulent  flow  with  chemistry. 
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Piguro  4,4,  Mach  Humber  Contours  for  SOME  notzle 
turbulent  flow  with  chemistry. 


